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SYMBOLS 


Sound  Speed 
Cycle  Number 

Slope  of  the  shock  velocity/particle  velocity  curve 
Time  Step 

Internal  Energy  Density 
Burn  Fraction 
Shear  Modulus 
Bulk  Modulus 
Pressure 

Hugonlot  Pressure 
Temperature  or  Time 
Shock  Velocity 
Particle  Velocity 
Yield  strength 
Distention  Ratio 

Grunelsen  ratio  or  ratio  of  specific  heats 

Grunelsen  ratio 

Volumetric  Strain,  1  -  po/p 

Compression,  p/po 

Excess  Compression,  p/po  -  1 

Poisson's  Ratio 
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TOODY  MODIFICATIONS 
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CONTROL  CARDS 

SOME  TIPS  ON  SUCCESSFUL  RUNNING 
REFERENCES 


APPENDIX  A 


SECTION  I 


INTRODUCTION 

^TOODY  Is  a  two  dimensional  (plane  or  axl symmetric) ,  multi- 
material,  Lagranglan  wave  propagation  computer  code.  It  was  written 
by  Sandla  Laboratories  In  the  late  1960'sj  References  1  and  2 
describe  one  of  the  early  versions.  fThe  Air  Force  Armament  Laboratory's 
(AFATL)  TOODY  Is  that  version  updated  extensively  by  both  Sandla  and 
the  Air  Forcej (primarily  Lt  Col  Osborn  while  at  AFATL  and  the  Air 
Force  Weapons  Laboratory).  • 

^The  purpose  of  this  report  Is  to  document  changes  made  to  the 
code  by  Sandla  and  the  Air  Force  and  to  provide  Input/ output  detail^/ 

It  Is  expected  that  the  reader  will  be  familiar  with  References  1  and 
2  prior  to  reading  this  report. 

AFATL  TOODY  runs  on  the  CDC  6600  computers  at  Eglln  Air  Force 
Base,  Florida  under  the  SCOPE  3.4  system.  It  uses  the  CDC  FTN  compiler 
and  the  CDC  UPDATE  system.  UPDATE  is  used  to  tailor  coding  to  the 
specific  problem  being  run.  It  runs  most  problems  at  a  speed  of 
1.5  x  10"3  CP  seconds  per  zone  per  cycle. 

The  code  can  be  made  available  by  contacting  Lt  Col  Osborn  at 
AFATL/DLJW,  Eglln  AFB,  FL  32542  (AV  872-2141). 

A  rezoner  exists  for  Sandla's  TOODY.  It  Is  planned  to  convert 
It  for  use  with  AFATL  TOODY.  The  conversion  has  not  yet  been  made. 

The  changes  required  Involve  provision  for  more  equations  of  state  and 
more  zone  variables. 


SECTION  II 
TOODY  MODIFICATIONS 

This  section  provides  Information  on  code  modification  made 


both  by  the  Air  Force  and  by  Sandla  Laboratories.  The  changing 
office  Is  Indicated.  The  modifications  are  not  discussed  In  any 
particular  order. 

EQUATIONS  OF  STATE 
ELASTIC/PLASTIC  (STATE  D-AIR  FORCE 
EOS  1  has  been  slightly  modified.  An  option  Is  now ‘available 
for  Hugonlot  pressure  as  a  function  of  u(excess  compression)  *  p/pQ-l . 

PH  =  poco2(y  +  kly2  +  •“  +  k5P6) 

Zone  variable  BT2  is  used  to  store  the  number  of  cycles  a  zone 

has  been  In  the  plastic  state.  Zone  variable  BT1  is  used  to  store 
the  number  of  cycles  a  zone's  pressure  has  been  less  than  the  input 
value  PMIN.  These  zone  variables  can  be  used  as  flags  for  Indicating 
plastic  and  spalled  zones  on  grid  plots. 

HIGH  EXPLOSIVE  (STATE  2)-AIR  FORCE 
An  option  has  been  added  so  that  the  Jones  -  Wilkins  -  Lee  (JWL) 
equation  for  detonation  products  can  be  used.  This  equation,  described 
fully  In  Reference  3,  Is  as  follows. 

P  =  A(1  -  )e-Rl/n  +  b(1  -  -^-)e"R2/ri  +  wpE 

where  A,  B,  R-| ,  R2  and  w  are  material  constants.  At  high  expansions, 
n  approaches  zero  and  P  approaches  wpE,  a  gamma-law  formulation. 
Coefficients  for  several  high  explosives  are  given  in  Table  I I -1 . 


JWL  COEFFICIENTS 


If  a  zone  has  not  burned,  Its  pressure  Is  computed  from 


P  =  K0y 
or 

P  =  Poco2n/0  -  sn)2 

In  the  latter  case  c0  and  s  are  set  to  values  appropriate  for 
Composition  B  —  c0  ■  2.88E5cm/sec  and  s  =  1.6.  In  either  event, 
a  message  Is  printed  warning  of  a  possible  low  order  detonation  If 
P  exceeds  25  k'ilobars  prior  to  detonation  time. 

Burn  time  options  have  been  increased  by  allowing  the  user  to 
set  a  constant  value  to  be  added  to  burn  time  computed  by  geometric 
means.  This  allows  more  accurate  burning  around  corners. 

IDEAL  GAS "SAND! A 

EOS  3,  the  Ideal  gas  equation  of  state  is  identical  to  Sandla's 
except  that  energy  is  now  automatically  added  on  cycle  0  to  cause  the 
Input  sound  speed  to  be  valid. 

COMPRESSIBLE  EARTH  MEDIA-AIR  FORCE 

EOS  4  has  been  written  to  provide  a  reasonable  formulation  for 
soil  and  porous  rocks.  In  addition,  a  built-in  concrete  model  Is 
available. 

EOS  4  assumes  that  pressure  Is  a  function  of  y  and  (the 
maximum  value  of  y  seen  by  the  zone).  The  specific  form  is  shown  in 
Figure  II-l.  An  energy  variation  on  pressure  is  not  considered. 

This  should  be  valid  for  all  but  very  high  energy  deposition  problems 
since  values  for  Grunelsen's  ratio  for  earth  media  and  concrete  are 
around  0.1  to  0.3  (as  opposed  to  values  from  1.5  to  3  for  metals). 


To  define  pressure  the  user  Inputs  pQ,  cp,  PMIN,  PELAS,  PMID,  yMID, 
PLOCK,  yLOCK ,  coL  and  s.  p0  Is  the  Initial  density  of  the  material. 
Cp  Is  Its  seismic  P-wave  velocity  (low  pressure  longitudinal  wave 
velocity).  PMIN  is  the  minimum  pressure  the  material  can  withstand 
without  breaking  up  in  tension.  As  seen  in  Figure  II-l,  Initial 
loading  is  along  a  curve 


P  =  poco  y 


Unloading  after  reaching  PELAS  but  prior  to  reaching  PLOCK  is 
accomplished  along  a  line 


P  -  Pymax  "  *^ymax  ^ 


where  P  is  the  pressure  at  u  =  umax  on  the  crushup  curve  and  K  Is 
Umax 

O 

a  slope  which  varies  from  p0cQ  when  y  =  yELAS  to  KLgCK  when  y  -  Mlqck 


Subsequent  reloading  is  also  along  this  line  until  the  crushup  curve 
Is  reached.  When  ymax  has  reached  ylOCK  all  loading,  unloading  and 
reloading  is  along  the  lockup  curve.  The  lockup  curve  is  described  by 
three  connecting  segments 

y  <  uL0CK’  P  =  PL0CK  '  KL0CK  ^LOCK 


u  >  wliOCK  ®UT  y  <  yHI’ 

P  s  poL  coLnL  I  0  ■  s\)2 

where  hl  =  1  -  p0l/p-  Thl s  is  a  form  suitable  for  describing  the 
Hugonlot  when  the  shock  velocity/partical  velocity  curve  of  the 
Hugonlot  can  be  described  by 

us  =  coL  +  sUp 

The  density  poL  is  computed  by  TOODY  to  insure  that  this  lockup  curve 
intersects  the  crushup  curve  at  (PLOCK,  yLOCK) . 


if  y  >  phi*  p  =  phi  +  ’Sii  -  yHr 

where  y^j,  PHI  and  K^j  are  calculated  by  TOODY  simply  to  prevent 
running  Into  the  pole  In  the  Hugonlot  form. 

Unloading  Is  limited  at  all  times  by  PMIN. 

The  yield  strength  of  the  material  Is  a  function  of  pressure. 

The  specific  form  provided  Is  shown  In  Figure  II-l .  If  P  <  PMIN,  Y  =  0 
If  P  >  PMIN  but  <  PV;  Y  is  calculated  by  interpolating  along  the  lines 
shown.  If  P  >  PV,  Y  is  set  to  Y V.  This  provides  a  Mohr-Coulomb 
curve  with  a  von-MIses  limit  of  YV.  Plasticity  is  checked  by  comparing 
Y  to  JIjJ.  If  exceeds  Y  all  deviators  are  multiplied  byjv/  3J^ 
This  is  the  same  flow  rule  used  for  metals  and  ignores  any  dllatancy 
In  the  material . 

Initial  bulk  sound  speed.  c0  is  computed  by  assuming 

c  2  =  (1  +  v)  c  2 

0  3TT  V)  p 

The  Rigidity  Modulus,  or  shear  modulus,  G,  Is  computed  from 

G  frr+v) K 

where  K  Is  the  current  bulk  modulus.  Provision  Is  made  for  limiting 
G  to  an  Input  value  G,^. 

Zone  variables  BT1  and  BT2  are  used  to  store  the  number  of  cycles 
that  P  Is  less  than  PMIN  and  number  of  cycles  during  which  the  zone 
Is  plastic. 

The  built-in  concrete  equation  of  state  is  described  in  Appendix  A 


to  this  report.  It  is  designed  for  5000  PSI  concrete  but  can  be  changed 
to  any  strength  by  changing  the  value  of  FPC  in  subroutine  SETCONC. 


Density  and  other  parameters  can  also  be  changed  in  the  subroutine. 

The  5000  PSI  model  has  been  used  to  compute  penetration  of  concrete 
at  penetrator  velocities  from  300  to  1600  fps.  It  works  very  well 
for  at  least  these  cases. 

A  model  has  also  been  developed  for  40  percent  porosity  dry  sand 
from  References  4  and  5.  The  Ru  relationship  is  shown  In  Figure  I I -2. 
The  sand  is  assumed  to  be  a  fluid  with  no  yield  strength.  Input 
parameters  (in  C6S)  are 

P0  1-6 
cp  0.6096E5 

pm  in  n 

v  0.5 
PELAS  50. E6 
PMID  1.E9 
pMID  0.2 
PL0CK  17.5E9 
uLOCK  0.425 
c6L  0.441E5 
S  2.33 

The  model  has  been  used  to  calculate  cavity  size  in  underground  HE 
detonations  and  has  provided  very  good  agreement  with  experiments  in 
the  Eglin  AFB  area. 

FOAM  (STATE  6>SANDIA 

EOS  6  provides  a  metal  foam  equation  of  state  routine.  It  is  the 


Sandia  ALPHA-E  model.  Ref  2  describes  the  ALPHA -P  model.  The  only 


difference  Is  that  In  the  ALPHA-E  model,  Pft  and  P$  can  be  functions  of 
the  Internal  energy  E.  The  specific  forms  provided  are 

pe  “  H<E>  peo 
ps  *  H(E)  Pe. 

where  H(E>  «  1  -  (k  +2)  (e-q)  +  (k+  1)  (|o)2,  K.  EQ,  Peo,  Ps0  are  all 
Input  quantities.  Generally,  k  should  lie  In  the  range 

2  <  k  <  -  1 


There  Is  no  strength  provided  In  this  model. 


GRAY  METAL  (STATE  7)-AIR  FORCE 


Equation  of  state  7  has  been  added  to  use  the  GRover  And  Young 
(GRAY)  metal  equations  of  state  developed  at  the  Lawrence  Livermore 
Laboratory  (Ref.  6).  GRAY  Is  a  three-phase  equation  of  state.  It 
is  coupled  with  State  1  rigidity  modulus  and  yield  strenth  options. 
Since  the  EOS  is  adequately  explained  In  Reference  6,  It  will  not 
be  discussed  here. 

Table  I I -2  lists  some  Input  values  required  by  the  EOS  and 
values  as  taken  from  Reference  6  for  several  metals.  Units  In  the 
table  are  as  they  should  be  for  Input  to  TOODY.  A  Mle-Grunelsen 
EOS  Is  used  for  the  solid  and  liquid  phases.  The  Input  variables 
c0  and  s  define  the  Hugonlot  and  yQ  and  a  define  the  Grunelsen 
parameter, 

y(p)  *  y0  -  an 

The  variable  ge  Is  the  electronic  energy  coefficient  In  Mbar  -cc/ 
mole-deg2.  The  melting  temperature,  Tm,  is  calculated  from  a 
Lindemann  law.  Tmo,  the  melt  temperature  in  deg  K  at  p  =  po,  is 
required  as  Input.  It  Is  shown  in  Reference  6  that  Tmo  is  approxi¬ 
mately  1.3  x  Tf^  where  T^  Is  the  normal  meltlna  temperature.  The 
variable  AW  is  the  atomic  weight.  Rj  defines  the  ratio  of  volumes 
Vj/Vq  where  Vj  is  the  volume  at  which  the  EOS  is  joined  to  a 
modified  van  der  Waals  EOS  for  the  vapor  region.  Rg  Is  Vg/V0  where 
VB  Is  the  excluded  volume  for  the  vapor  phase.  The  variable  AY  Is  a 
van  der  Waal's  coefficient.  TH  Is  a  join  parameter  explained  in  Ref  6. 


In  State  7,  GRAY  has  been  joined  with  the  rigidity  modulus  and 
yield  strength  options  of  State  1.  If  an  Input  variable  TRFLAG 
Is  set  to  1,  then  both  Y  and  PMIN  are  reduced  linearly  to  0  at 
melt  temperature. 

State  7  uses  the  BT1  zone  variables  to  Indicate  the  phase  of 
the  material.  BT1  =  0  means  solid,  1  means  melting,  2  means  liquid 
and  3  means  vapor.  BT2  contains  the  ratio  of  current  temperature 
to  melt  temperature. 


SLURRY  EQUATION-AIR  FORCE 


Equation  of  State  9  provides  a  detonation  products  equation 
suitable  for  slurry  explosives.  It  was  developed  by  personnel  at 
the  Lawrence  Livermore  Laboratory.  It  employs  two  JWL  equations. 
About  25%  of  the  slurry's  total  chemical  energy  is  released  upon 
detonation  into  the  first  JWL  equation.  Energy  is  then  released 
on  a  time  ramp  into  both  JWL  equations.  During  this  phase  the 
pressure  contributed  by  each  equation  is  proportional  to  the  amount 
of  energy  released  so  that  when  all  of  the  energy  has  been  released 
pressure  is  provided  by  the  second  JWL  equation  only.  The  equations 
for  State  9  are 

P  =  F  (1  -  FR)Pi  +  FR  •  P2 

Where  F  is  the  burn  fraction  used  in  the  first  energy  release  and 
FR  Is  a  burn  fraction  used  to  release  energy  along  a  time  ramp. 

P-|  is  the  pressure  from  the  first  JWL  equation  and  P2  is  the  pressure 
from  the  second  JWL  equation.  These  pressures  are  calculated  from 

Pi  =  a1  (1  -lJ#-)e  "Rli/fn  +  B1  (1  "R2i/n  +  wiPE 

Where  1  =  1  or  2  and  n  =  p/pQ 

During  the  ramp  burn,  internal  energy  density  E  is  increased  by 
EZ  ■  EZDEL  •  DT/TRAHp 

where  EZDEL  is  the  amount  of  energy  to  be  deposited  on  the  ramp,  DT 
Is  the  time  step  and  T^p  is  the  time  spent,  EZ  is  not  allowed  to 
exceed  EZDEL.  The  time  step  is  restricted  during  ramp  burn  to  be  no 
larger  than  0.05DT^mp- 


The  first  burn  fraction,  F,  Is  computed  as  for  any  explosive. 
FR  Is  computed  as  follows. 

FR  «  0  If  F  <  1 
FR  =  MX  (1,  EZ/EZDEL)  If  F  -  1 
State  9  uses  zone  variables  BT1  and  BT2  as  follows.  BT1  Is  FR. 

BT2  Is  the  amount  of  energy  released  to  the  current  time,  EZ. 


ZONE  VARIABL£S“*IR  FORCE 

The  zone  variables  have  been  changed  somewhat  from  References 
1  and  2.  Three  variables,  BTT  ;  BT2  and  BT3  have  been  added  bringing 
the  total  to  22.  The  deviatorlc  viscosity  option  has  been  dropped 
and  those  variables  DQXX,  DQXZ  and  DQZZ  are  used  to  store  strain 
values.  The  zone  variable  lineup  now  Is 

MAT  Material  indicator 

X  Value  of  X  coordinate 

Z  Value  of  Z  coordinate 

UX  Velocity  In  X  direction 

UZ  Velocity  In  Z  direction 

TXX  Stress  components.  TYY 

TYY  Is  hoop  stress 

TZZ 
TXZ 

Q  Artificial  viscosity 

DQXX  Strain  Component.  EXX 

DQZZ  Strain  Component,  EZ2 

DQXZ  Strain  Component,  EXZ 

RHO  Density 

C  Longitudinal  Sound  Speed 

E  Internal  energy  density 

AL  Area  of  zone  in  X  -  Z  plane 


16 


XM 


Mass  of  zone 


BT  Burn  time  If  IES  ■  2 

Plastic  work/unit  volume  If  IES  *  1 
and  work  hardening  is  being  used. 

If  IES  *  1  and  work  hardening  not  used, 
then  BT  contains  energy/unit  mass  to  be 
deposited  In  the  zone  In  time  TDEP. 

If  IES  »  6,  BT  contains  value  of  ALPHA, 
the  distention  ration 

BT1  Number  of  cycles  a  zone's  pressure  has 
had  to  be  reset  to  PMIN 
If  IES  *  1  or  4 

BT2  Number  of  cycles  a  zone  has  been  plastic 
IF  IES  =  1  or  4 

BT3  Hoop  strain  component,  EYY 
For  boundaries  (MAT  >  20),  only  MAT,  X,  Z,  UX  and  UZ  are  non¬ 
zero  with  the  exceptions  that  for  MAT  *  22  the  stress  components 
carry  the  applied  stresses  and  for  boundaries  along  IT  and  JT  only 
MAT  and  XM  are  non-zero  (XM  Is  -1  for  these  zones). 

Zone  positions  and  velocities  are  valid  at  the  II,  JJ  Intersection 


Other  zone  quantities  apply  to  the  zone  contained  between  II-l  and  II 
and  JJ-1  and  JJ. 


Strain  components  are  not  needed  in  TOODY  to  calculate  stresses. 
They  are  saved  merely  as  a  convenience  for  plotting  and  editing 
purposes . 

RESTARTING  -  AIR  FORCE 
RESTART  FILE  FORMAT 

The  restart  file  format  has  been  changed.  It  is  as  follows 
RECORD  VARIABLES 

1  ID ( 1 )  through  ID(9),  NCYCLE ,  T,  DT,  IACT.  IT. 
IHULLIN.  ISPEC,  NOVAR,  IJ(1)  through  IJ(lO),  JMIN(l) 
through  JMIN(IO) ,  JMAX(l)  through  JMAX(IO). 

2  JACT,  JKACT  (1) 

If  JACT  GT  0  then  the  next  record  appears  on  the  tape 

3  JKACT(2)  through  JKACT(IT) 

4  NV,  NOMAT,  IALPH 

If  IHULLIN  GT  0  then  the  next  four  records  appear  on  the  tape 

5  RP0S(1)  through  RP0S(200),  YP0S(1)  through  YP0S(200) 
NSTN,  NREC 

6  NSTPT (11)  through  NSTPT(400,4) ,  KBNDRY ,  KBDY(1,1) 
through  KBDY(IT  JT)  XAA(1  1)  through  XAA(400,  4), 
ZAA( 1 ,  1)  through  ZAA(400.  4) 

7  HPEAK(1  1)  through  HPEAK(200,4\  HCURR(1,  1)  through 

HCURR ( 200  4) 

PPRESS ( 1  1.1)  through  PPRESS(200  2,4).  TPRESS(2), 
NPRESS1  ,  NPRESS2  THULL 
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If  ISPEC  GT  0  then  the  next  two  records  appear  on  the  tape 

1  ITREC.PTLAST(IOO)  TLAST,  EXX(IOO).  EXZ(IOO) 

EZZ(IOO) ,  ITNAT,  NOSTA,  ITPOS,  MPTS(IOO).  XXTO(IOO) 
ZZTO(IOO) ,  EYY(IOO) 

2  IITIME(1 )  through  IITIME(ISPEC), 

JJTIME(l)  through  JJTIME( ISPEC), 

XTIME(l)  through  XTIHE  (ISPEC), 

ZTIME(l)  through  ZTIME  (ISPEC), 

INUM(l)  through  INUM  (10) 

Each  II  column  Is  a  record.  II  =  1  through  II  =  IT  follow.  Each 
record  contains  the  22  zone  variables  for  JJ  *  JMIN  (JTI)  through 
JMAX(JTI)  where  JTI  Is  the  section  number  for  the  II  column. 

The  last  record  on  the  dump  Is  an  end-of-file  Each  dump  therefore 
Is  separated  by  an  end-of-file, 

RESTARTING  PROCEDURES  -  AIR  FORCE 

If  RSMIN  Is  greater  than  zero,  TOODY  will  create  restart  dumps 
based  on  elapsed  CP  time  or  on  problem  time.  Each  dump  Is  put  on 
TAPE  25  In  the  Initial  run  of  the  problem.  Up  to  20  files  are 
allowed  on  TAPE  25.  Subsequent  dumps  are  put  on  TAPE  26  and  then 
TAPE  27.  After  20  dumps  have  been  put  on  TAPE  27,  TAPE  25  will  be 
called  again.  There  Is  no  built-in  tape  length  counter.  For  very 
large  problems  the  20  dumps  allowed  may  be  too  many.  If  so:  this 
must  be  changed  at  UPDATE  line  AFATL3.1  in  the  main  routine.  Normally, 
dumps  are  required  so  seldom  that  TAPE  25  is  never  filled. 


When  restarting,  TAPE  25  must  be  positioned  to  the  beginning 
of  the  correct  dump.  This  can  be  accomplished  In  many  ways  -  e.g., 
by  using  the  correct  number  of  COPYBF's  in  the  control  card  deck. 

After  restarting,  further  dumps  are  put  on  TAPE  26  (up  to  20)  and 
then  TAPE  27.  Finally,  TAPE  25  will  be  called  again. 

To  restart,  the  user  puts  in  a  blank  ID  line  as  the  first  card 
In  Input.  Other  cards  should  be  Input  also.  Changes  can  be  made 
in  non-zone  variables  such  as  material  properties,  TMAX,  plotting 
input,  etc.  Zone  variables  such  as  material  property  indicators 
and  mesh  sizes  can  be  changed  only  for  columns  greater  than  IACT  +  1 
•  ZONING  OPTIONS -AIR  FORCE  AND  SANDIA 

Input  cards  for  each  geometry  region  are  as  follows: 

ISET 

XO  XC0N(1),  XCON ( 2 ) ,  XC0N(3).  XC0N(4).  JMIN.  JMAX 
ZO  ZC0N(1),  ZC0N( 2 ) ,  ZC0N(3).  ZC0N(4).  IMIN.  IMAX 
ISET  =  2,  3,  4,  5,  and  6  have  not  been  changed  from  References  1  and 
2.  ISET  =  1  has  been  changed  somewhat.  In  this  option 

X  =  XO  +  XCON ( 1 ) ( JJ  -  JMIN)  +  ( XCON ( 2 )  +  XCONf 3) ( JJ  -  JMIN) 

+  XCON (4)  (Z  -  ZO))  (Z  -  ZO)  +  XCON  (5)  (II  -  IMIN) 

Z  *  EZ  +  (ZCON  (1)  +  ZCON ( 2 )  (JJ  -  JMIN))  (II  -  IMIN)  +  ZC0N(5)(JJ  -  JMIN) 

where  EZ  =  ZO  if  ZC0N(3)  =  0 

and  EZ  =  ZO  +  0.5  ZCON ( 3 ) ( I I  -  IMIN)  (II  -  IMIN  -  1) 

if  ZCON (3)  t  0 

In  the  ISET  =  1  option,  XC0N(5)  is  read  in  the  JMAX  field  and  ZC0N(5) 
is  read  in  the  IMAX  field. 


Sandla  has  added  an  ISET  =  7  option.  In  this  option  the  bound¬ 
aries  of  the  region  between  IMIN,  JMIN,  IMAX  ,  JMAX  Indices  can  be 
either  straight  line  segments  or  arcs  of  circles,  with  the  restrict¬ 
ion  that  no  arc  may  include  an  angle  of  180°  or  larger.  Mesh  nodes 
on  a  line  of  constant  JJ  will  be  such  that  the  arc  length  from  the 
point  at  IMIN  is 
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where  SIT  is  the  total  arc  length  of  the  line  of  constant  JJ. 
Similarly  nodes  on  a  line  of  constant  II  will  be  such  that 
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The  boundary  arcs  may  be  concave  or  convex.  If  they  are  concave,  the 
input  value  of  the  radius  is  negative.  For  straight  line  boundaries, 
the  Input  Is  a  0  radius.  Three  Input  cards  are  required.  Input  is 
Xq,  X!*,  X2*,  X3*,  X4*,  JMIN,  JMAX 


Z0,  Z-j*,  Z2*,  Z3*,  Z4*,  IMIN,  IMAX 
R12’  R23’  r34>  r41  »  eXI>  eXj 

If  eXI  or  eXJ  are  zero  then  the  exponent  will  be  a  1.  See  Figure  I I -3 . 
ISET  =  8  is  not  used 

ISET  *  9  and  10  options  have  been  added  by  the  Air  Force.  ISET  =  9 
is  used  for  regions  with  circles  or  Z  *  constant  lines  along  II  columns 
and  X  =  constant  lines  along  JJ  rows. 

Input  is  as  follows: 

X0  Is  X  for  center  of  II  -  IMIN  circle 
XC0N(1)  is  Z  for  center  of  II  =  IMIN  circle 
XC0N(2)  is  radius  of  the  II  =  IMIN  circle 
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If  II  =  IMIN  is  a  straight  line  then  the  constant  Z  value  for 
this  line  Is  XC0N(1).  In  this  case  XC0N(2)  must  be  zero. 

XC0N(3)  is  the  value  of  X  for  JJ  =  JMIN 

XC0N(4)  Is  the  value  of  X  for  JJ  a  JMAX 

ZO  is  X  for  center  of  II  -  IMAX  circle 
ZC0N(1)  Is  Z  for  center  of  II  =  IMAX  circle 
ZC0N(2)  is  radius  for  the  circle  at  II  =  IMAX 

If  II  =  IMAX  is  a  straight  line,  then  ZCON ( 2 )  must  be  zero  and  II*IMAX 

Is  a  straight  line  along  Z  =  ZC0N(1). 

ZC0N(3)  is  any  negative  number  If  the  center 

of  the  IMIN  circle  Is  left  of  th*  II  *  IMIN 
circular  surface. 

A  positive  or  zero  value  Indicates  the 
center  is  to  the  right  of  the  II  =  IMIN 
surface. 

ZCGN(4)  Negative  if  center  of  IMAX  circle  is  left 
of  the  II  *  IMAX  surface. 

Positive  or  zero  if  center  is  to  the  right. 

ISET  =  10  Is  an  additional  shape  option  added  to  allow  elliptical 
zoning.  The  ellipse  must  be  centered  at  (0,  0)  and  a  full  half  of 
the  ellipse  must  be  modelled.  The  inner  surface  is  II  =  IMIN  and 
has  axes  XC0N(1)  in  the  X  direction  and  XC0N(2)  in  the  Z  direction. 

The  outer  surface  Is  II  =  IMAX  and  has  axes  XCON { 3 )  and  XC0N(4)  In 
the  X  and  Z  directions  respectively.  If  ZC0N(1)  is  non-zero  then 


II  =  IMAX  is  drawn  with  axes  equal  to  those  IMIN  +  ZC0N(1)  —  l.e., 
ZC0N(1)  is  the  thickness  In  this  case. 

ISET  =  9  and  ISET  =  10  are  also  illustrated  in  Figure  II  -  3. 
Initial  velocity  values  are  no  longer  input.  Zero  values  are 
assumed  for  all  material*.  If  non-zero  values  are  desired  they  must 
be  UPDATED  into  the  SETUP  1  subroutine. 

APPLIED  PRESSURE  BOUNDARY-AIR  FORCE 
The  original  Sandia  TOODY  offered  a  time  varying  applied  pressure 
form.  This  has  been  taken  out.  If  applied  pressure  is  desired  it 
can  be  taken  from  a  HULL  time  history  tape  or  a  subroutine  must  be 
UPDATED  into  the  code.  The  subroutine  must  begin  with 
SUBROUTINE  PRESSIN  (TIMEN, IPARG, PRESS) 

TIMEN  and  IPARG  are  transferred  Into  the  routine  and  PRESS  must  be 
computed  within  the  routine.  TIMEN  is  the  problem  time.  IPARG  Is  the 
zone  number  -  1,  2,  3,  or  4  (see  Reference  1).  PRESS  must  be  negative 
If  compression  is  desired  and  positive  for  tension.  This  is  true 
because  PRESS  Is  really  used  as  a  stress  input. 

The  largest  change  in  applied  pressure  is  that  used  to  take 
stress  values  from  a  HULL  time  history  tape  (TAPE  9)  and  apply  them 
to  a  TOODY  problem.  Also,  since  TOODY  and  HULL  time  history  tapes 
are  Identical  in  format  this  means  that  a  TOODY  problem  could  be 
broken  into  separate  problems  connected  by  a  pressure  surface.  It  Is 
Implied  in  using  this  routine  that  shock  impedance  differences  and 
TOODY  boundary  zone  distortions  are  such  that  important  physics 
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interactions  are  not  lost  in  the  decoupling. 

The  changes  in  TOODY  required  to  use  a  HULL  tape  as  input  are 
quite  extensive.  The  changes  are  not  as  general  as  would  be  desirable. 

A  small  amount  of  UPDATING  must  be  employed  to  use  HULL  Input. 

The  HULL  time  history  tape  format  is  as  follows. 

RECORD  VARIABLES 

1  Two  tape  numbers  (Ignored  by  TOODY) 

2  HULL  ZBLK  data  (200  words).  TOODY  uses 
NOP  =  the  number  of  particles  In  the  HULL 

problem  (particles  can  be  tracers  or 
time  history  station  locations) 

NPP  =  the  number  of  position  values  needed 
to  describe  a  particle's  location 
(2  for  two  dimensions) 

NSTN  *  The  number  of  time  history  stations. 

3  This  record  contains  (R,  Y)  position  pairs 
for  the  NOP  particles  and  Is  NOP  x  NPP  long. 
NOTE:  TOODY  assumes  NOP  x  NPP  is  no  larger 
than  1024.  If  more  room  Is  needed,  then 
the  dimensions  on  the  array  AHULL  must  be 
increased  in  the  HULLII  common  block.  The 
dimension  in  AHULL  must  never  be  less  than  1024. 
Lower  bits  of  the  R  and  Y  words  are  examined 

to  determine  if  the  particle  Is  a  station  and, 


if  so,  Its  number.  TOODY  also  assumes  no 
more  than  200  stations.  If  more  are  needed 
RPOS  and  YPOS  dimensions  must  be  changed  In 
HULLII. 

4  Each  record  from  record  number  4  on  contains 

time  history  data.  The  first  word  Is  T. 

(time)  with  NW  (the  number  of  words  of  data 
for  that  time)  masked  into  the  lower  bits  of 
T.  Since  there  are  15  variables  saved  per 
station,  the  number  o*  stations  with  data 
at  time  T  is  NW/15.  Masked  into  the  lower 
bits  of  the  first  of  the  15  variables  is  the 
station  number.  Stresses  are  contained  in 
variables  7  (radial  stress),  8  (Shear  stress), 
9  (axial  stress)  and  10  (hoop  stress). 

Subroutine  POSREAD  Is  called  from  PRESET  and  reads  the  first  3 
HULL  records  if  TOODY  Input  variable  IHULLIN  is  Set  to  1.  It  adds 
RHOFF  and  YHOFF  to  HULL  R  and  Y  station  positions.  RHOFF  and  YHOFF 
are  input  to  TOODY  and  are  offsets  between  the  HULL  and  TOODY  Coordin¬ 
ate  systems. 

Subroutine  PREAD  fills  the  PPRESS  (200,  2,  4)  and  TPRESS(2) 
arrays  from  HULL  data.  PPRESS  contains  station  pressures  at  two  times, 
TPRESS(l)  and  TPRESS(2).  These  two  times  are  such  that  TOODY  time  is 
contained  at  a  point  between  them.  The  first  dimension  in  PPRESS  is 
the  HULL  station  number.  The  second  identifies  the  time  -  either 


TPRESS(l)  or  TPRESS(2).  Variables  NPRESS1  and  NPRESS2  are  used  In 
a  flip-flop  manner  to  keep  up  with  time.  The  earliest  time  Is  contained 
In  TPRESS(NPRESSl).  The  latest  time  Is  In  TPRESS(NPRE$S2).  Initially 
NPRESS1  *  1  and  NPRESS2  -2.  To  avoid  copying  of  PPRESS  from  one  time 
to  the  next,  NPRESS1  and  NPRESS2  are  flip-flopped.  When  a  new  time  Is 
needed  It  and  stress  data  are  read  Into  TPRESS  and  PPRESS  at  NPRESS1 > 
since  that  data  Is  no  longer  needed.  Then  NPRESS1  and  NPRESS2 
exchange  values.  The  second  dimension  In  PPRESS  Identifies  time  and 
Is  either  NPRESS1  or  NPRESS2.  The  third  dimension  Identifies  the 
stress  component.  Radial  stress  Is  1,  shear  stress  Is  2,  axial  stress 
Is  3  and  hoop  stress  Is  4. 

T000Y  must  Identify  the  HULL  stations  (one  for  corners  and  two 
otherwise)  which  define  pressure  for  each  applied  pressure  zone. 

This  coding  has  not  been  generalized  and  must  be  rewritten  for  each 
problem.  Each  TOODY  zone  requiring  pressure  Input  Is  identified  In  an 
array  KBDY(II,  JJ).  For  any  (II,  JO)  the  value  of  this  array 
KBNDRY  =  KBDY  (II,  JJ)  Is  a  key  number  needed  to  Identify  the  HULL 
stations  providing  input  to  that  TOODY  zone.  KBDY  Is  dimensioned 
to  (IT.  JT)  by  the  preprocessor.  The  II,  JO  values  In  KBDY  must  be 
those  of  the  zone  being  processed  through  the  momentum  routine  In 
TOODY  MAIN.  These  are  zones  which  define  the  physical  boundary  in  the 
applied  pressure  region.  The  HULL  stations  Identified  with  this 
boundary  point  are  contained  In  the  array  NSTPT  (400,  4).  Two 
stations  are  encoded  In  a  215  format  In  the  NSTPT  word.  The  limit 


of  400  means  there  can  be  no  more  than  400  TOODY  zones  requiring  HULL 
input.  The  second  dimension  in  NSTPT  is  IPAR6  and  identifies  the 
surrounding  zone  as  1,  2,  3  or  4  (See  Reference  1 ) . KBNDRY  Is  also 
used  to  find  the  average  X  and  Z  values  in  TOODY  along  the  boundary 
of  the  zone  subjected  to  applied  pressure.  These  values  are  contained 
in  XAA  (400,  4)  and  ZAA  (400,  4).  The  second  dimension  is  again  IPAR6. 
These  numbers  are  required  because  TOODY  weights  input  from  the  two 
HULL  stations  by  using  range  from  (XAA,  ZAA)  to  the  HULL  station 
positions.  Initial  values  of  XAA  and  ZAA  are  saved  at  setup  time  and 
used  for  weighting  so  that  distortions  of  the  TOODY  grid  cannot  change 
the  weighting  functions. 

It  is  necessary  to  consider  an  example  to  understand  how  this  fairly 
complex  system  really  works.  Suppose  TOODY  is  modelling  a  cylinder 
subject  to  pressure  from  a  HULL  calculation.  The  mesh  for  this  is 
shown  in  Figure  II-4a.  It  is  7cm  long  and  5cm  in  radius.  TOODY  models 
just  one  half  of  the  cross-section  of  the  cylinder  and  sets  IALPH  =  2 
so  that  the  Z  axis  is  an  axis  of  rotation.  The  physical  boundaries  of 
the  cylinder  are  at  II  =  1,  II  =  6  and  JJ  =  1  and  JJ  =  4.  Material 
numbers  are  shown  in  II  -  4b.  The  material  indicator  input  cards  for 
this  problem  would  be 


Symmetry  (MAT  =  23)  is  used  alone  the  centerline  ( JJ  =  l) .  The  other 
boundary  zones  are  applied  pressure  (MAT  *  22).  IT  is  7  and  JT  Is  5 
because  one  additional  II  and  one  additional  JJ  line  must  be  provided 
to  carry  boundary  indicators.  The  II  =  6  and  JJ  =  4  physical  boundary 
lines  contain  the  material  number  for  the  actual  cylinder  material. 

The  momentum  routine  will  process  11=1  throuqh  11=6.  For  each  II 
value  it  will  process  JJ  values  from  1  through  4.  It  will  not  look 
at  II  =  7  or  JJ  =  5  lines.  For  each  (II.  JJ)  the  momentum  routine 
requires  values  of  stress,  density,  etc  for  the  four  zones  surrounding 
the  point.  This  is  shown  in  Figure  II  -  4c.  The  zones  surrounding 
(II,  JJ)  are  numbered  as  shown  in  the  figure.  The  dashed  line  indicates 
the  boundaries  of  the  diamond-shaped  momentum  zone.  If  all  four 
zones  contain  real  material  -  in  this  case  material  1  -  then  the  applied 
pressure  routine  (ENTRY  PRESURE  in  SUBROUTINE  OUTPUT)  is  not  called. 

If  any  of  the  zones  contain  a  22  indicator  then  it  is  called  with 
IPARG  being  the  zone  number  of  the  zone  requiring  pressure  data.  For 
example,  consider  that  TOODY  is  processing  the  point  II  =  1,  JJ  =  3. 

Then  material  numbers  in  the  four  zones  are  as  shown  in  Figure  II  -  4d. 
Zones  1  and  4  contain  real  material.  Zones  2  and  3  are  boundary  zones 
requiring  applied  pressure.  PRESURE  will  be  called  with  IPARG  =  2 
and  then  it  will  be  called  again  with  IPARG  =  3.  Figure  II  -  4e 
shows  another  possible  situation.  In  this  case  II  =  1  and  JJ  =  4  -  a 
corner  point.  Data  is  required  in  this  case  for  three  zones,  IPARG  =  1, 
2  and  3.  Along  JJ  -  4  data  will  be  needed  for  zones  1  and  2  (See 
figure  II  -  4f).  At  the  other  outer  corner,  II  =  6  and  JJ  =  4, 


pressure  Is  needed  for  zones  1,  2  and  4  (Figure  II  -  4g).  Along 
II  ■  6,  at  mid  JJ  values  such  as  3,  data  is  needed  for  zones  1  and 
4  (Figure  II  -  4h).  Pressure  data  will  not  be  required  along  JJ  «  1 , 
the  bottom  boundary,  except  at  the  corners.  For  example,  at  II  *  1, 

JJ  s  1  material  indicators  are  as  shown  in  Figure  II  -  4i.  Since  II  =  1 
is  an  axis  of  symmetry,  PRESURE  will  be  called  for  both  zone  2 
(IPARG  =  2)  and  zone  3  (IPARG  =  3).  If  there  is  no  shear  stress,  then 
the  stress  transferred  is  the  same  for  both  zones.  The  shear  stress 
component  must  be  reversed  in  sign  if  it  is  non-zero.  The  same  comments 
apply  to  the  corner  at  II  =6  and  JJ  =  1  with  the  exception  that  the 
zones  have  material  indicators  as  shown  In  Figure  II  -  4j. 

We  have  identified  the  zones  in  this  sample  problem  requiring 
applied  pressure  input  from  HULL.  The  user  must  provide  coding  to 
identify  the  NSTPT,  XAA  and  ZAA  arrays.  SUBROUTINE  SETUP  contains 
a  routine  to  provide  these  needed  variables  for  a  particular  problem. 

It  can  be  modified  or  a  new  routine  written.  For  illustration  the 
coding  required  to  modify  it  for  this  sample  problem  will  be  discussed. 
HULLIN.55  through  HULLIN.84  can  be  deleted  and  replaced  by 

DO  5900  11=1,  ITEND 
ICT  *  NV  *  JMAX(l) 

DO  5003  IPP  =  1  ,  ICT 

5003  IC  (IPP,  1)  =  IC( IPP ,  2) 

DO  5004  IPP  =  1,  ICT 

5004  IC  (IPP,  2)  =  IC  (IPP,  3) 

BUFFER  IN  (NTP1 ,  1)(IC(I,3),  IC(ICT,3) 

IF  (UNIT  (MTP1 ) )  5005,  5005,  5005 


5005  CONTINUE 


This  coding  will  put  the  II  -  1  column  in  IC  (  -,  1),  II  in  IC  (  -,  2) 
and  II  +  1  in  IC  (  3)  for  each  II  value.  HULLIN.85  through  HULLIN.90 

can  be  replaced  by 

JT1  =  JMAX(l) 

DO  5890  00  =  1 ,  JT1 

<m  - 

M  =  NV  *  (JJ  -  1)  +  1  t, 

IF  (II  .  EQ.l )  GO  TO  5006 
IF  (JJ  .  EQ.4)  GO  TO  5006 

IF  (II  .  EQ.6)  GC  TO  5006  - 

GO  TO  5890 

For  each  (II,  JJ)  point  requiring  NTSPT,  XAA  and  ZAA  data,  the  coding 

V 

sends  the  computer  to  statement  5006.  Statement  5006  will  begin  ► 

special  coding  to  provide  NTSPT,  XAA  and  ZAA  values  for  IPARG  =1 
through  4  as  needed.  HULUN.  101  through  HULLIN. 132  must  be  replaced. 

(- 

This  coding,  inside  a  DO  LOOP  on  IPARG  provides  average  X  and  Z  values 
for  each  boundary  zone,  XAVG  and  ZAVG,  and  also  provides  X  and  Z  values, 

XXI  and  ZZ1 ,  which  lie  along  the  boundary  line.  XAVG  and  ZAVG  will  be 
stored  later  as  XA7  <nd  ZAA.  The  (XXI,  ZZ1)  values  will  be  used  to 
define  a  normal  to  the  boundary  line.  This  normal  is  needed  to  allow 
the  code  to  find  the  nearest  two  HULL  stations.  Subsequent  coding  will 

i 

assume  that  corners  (where  only  one  station  is  desi.'ed)  will  be 
identified  by  values  of  XXI  -  777.  The  replacement  coding  could  be 

GO  TO  (5010,  5020,  5030,  5040), IPARG 

•  > 

5010  IF  (II.  LE.5.  AND.  JJ.  E0.4)  5011,  5012 

5011  XAVG  =  0.5*( X(M,  2)  +  X(M,3)) 

ZAVG  -  0.5*  (Z(M,  ?)  +  Z(M,  3)) 

I 

XXI  =  X(P,  3) 

34 
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ZZ1  =  Z{M,  3) 

GO  TO  5300 

5012  IF  (II.  EQ.  6.  AND.  JJ.  EQ.  4)  5013,  5014 

5013  XAVG  =  X(M,2) 

ZAVG  =  Z(M,  2) 

XXI  =  777 

GO  TO  5300 

5014  IF  (II.  EQ.6)  5015,  5300 

5015  XAVG  =  0.5*  (X(M,2)  +  X(M  +  NV,  2)) 

ZAVG  =  0.5*  (Z(M,  2)  +  Z(M  +  NV,  2)) 

XXI  =  X(M  +  NV,  2) 

ZZ1  =  Z(M  +  NV,  2) 

Coding  for  IPARG  =2,3  and  4  would  follow  but  will  not  be  shown  here. 
The  coding  above,  combined  with  that  already  in  SETUP  will  send  the 
code  to  statement  5300  with  (XAVG,  ZAVG)  and  (XXI,  ZZ1)  values  set  to 
those  needed  for  each  physical  boundary  zone  and  each  value  of  IPARG. 
The  remaining  coding  in  SETUP  will  handle  the  rest  of  the  problem 
automatically. 

The  only  further  complication  occurs  with  slide  lines.  Slide 
lines  cause  two  problems.  Slide  lines  car.  be  the  last  line  In  a 
section,  II  =  IJ  (JTI),  or  the  first  line  In  the  next  section, 

II  =  IJ(JTI)  +  1  (the  latter  is  recommended).  Whichever  Is  used, 
the  slide  line  Is  co-located  In  space  with  another  line  which  provides 
the  other  side  of  the  slip  surface.  This  can  complicate  the  SETUP 
coding  needed  to  set  the  NTSPT,  XAA  and  ZAA  array  values.  The 
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second  problem  occurs  if  a  slide  line  or  the  other  side  of  the  slip 
surface  actually  provides  an  applied  pressure  surface.  For  example, 
consider  Figure  I I -5 .  In  this  case  the  slide  line  is  II  =  4  and  its 
neighbor  II  =  3  has  a  portion  of  itself  extending  into  the  HULL  flow 
field.  The  T000V  slide  routine  was  written  by  Sandia  to  be  very 
general,  but  does  not  properly  handle  this  case  without  modifications. 
When  processing  II  =  3  and  JJ  -  1  through  4  zones,  TOODY  may  assume 
a  free  surface  boundary  since  o,.'  side  of  the  slide  line  extends 
further  than  the  other.  This  can  be  fixed  quite  simply  by  adding  a  few 
lines  of  coding  to  the  SLIDE  subroutine.  In  this  case,  coding  can  be 
added  In  ENTRY  SLIDE  2.  The  coding  to  be  added  would  be 
*  DELETE  T00DY3. 5280 

IF  (DO(JJ)  -  D1MIN.LE.ZER0)  K4  =  IC(M,IP)  =  22 
and  *  DELETE  TOODY 3.  5232 

IF  (D0(JJ  +  1)  D1MIN.LE.ZER0)  K1  =  IC(I,IP)  =  22 
DO  and  D1MIN  are  discussed  in  Reference  1. 


PLOTTING-AIR  FORCE 


Grid,  velocity  vector  and  zone  variables  vs  X  type  plotting 
routines  have  been  built  into  TOODY.  The  routines  use  only  a  small 
amount  of  central  memory  and  are  highly  flexible.  In  addition,  they 
can  be  easily  changed  to  provide  specialized  plots  such  as  warhead 
fragment  angles  vs  radial  position  of  the  fragments.  They  are 
written  in  SC  4020  language. 

The  plots  are  divided  into  two  general  types.  The  input  variable 
I GRID  provides  GRID  plots  which  include 

(1)  Distorted  grid  plots 

(2)  Velocity  vector  plots  with  material  and  physical 
boundaries  shown 

(3)  Material  and  physical  boundary  plots  with  BT1  or  BT2 
values  printed  in  zones  in  which  these  values  exceed  0. 

The  routine  can  be  easily  modified  to  print  peak  pressure,  maximum 
tensile  strain,  etc.  The  second  type  of  plot  is  controlled  by  the 
input  variable  IGRAF  and  provides  GRAPH  plots.  These  are  plots  of 
zone  variables  vs  X.  One  plot  can  contain  plots  for  several  II 
columns.  The  plots  are  restricted  to  plots  vs  X  because  of  the  way 
in  which  data  is  stored  on  intermediate  files  TAPE  21  and  TAPE  22. 
Generally  X  increases  with  JJ.  The  plotter  brings  in  one  II  column 
at  a  time  and  can  easily  plot  any  variable  as  it  changes  with  JJ 
values.  To  plot  along  a  constant  JJ  row  with  II  varying  would 
require  an  inordinate  amount  of  computer  time.  These  routines  will 
now  be  discussed  in  detail  and  examples  shown. 


GRID  PLOTTING 


The  total  number  of  grid  plots  to  produce  at  any  time  Is  IGRID. 
Proglem  times  for  plotting  are  controlled  by  TGRDDEL.  It  Is  assumed 
that  plots  are  desired  at  T  s  0,  T  *  TGRDDEL,  T  B  2  TGRDDEL,  .... 
through  the  end  of  the  problem.  TGRDMIN  can  be  used  to  start  grids 
at  a  problem  time  other  than  0  and  TGRDMAX  can  be  used  to  cease 
plotting  at  a  certain  problem  time.  These  are  now  set  automatically 
at  0  and  1E6  seconds  respectively.  Changes  to  these  values  must  be 
made  by  UPDATING  new  values  Into  subroutine  READGRD  (which  reads  In 
GRID  plotting  control  variables).  Input  variables  for  each  plot 
desired,  from  1  to  IGRID,  are  on  two  cards.  The  first  Is 

IGRDMIN,  IGRDMAX,  JGRDMIN,  JGRDMAX ,  I0PTI0N 
the  second  card  is 

XGRDMIN,  XGRDMAX,  ZGRDMIN,  ZGRDMAX ,  XHORIZ 
IGRDMIN,  IGRDMAX,  OGRDMIN  and  JGRDMAX  are  bounding  values  of  II  and  JJ 
to  be  included  In  the  plots.  These  variables  allow  the  user  to  plot 
specific  portions  of  the  grid  —  he  need  not  always  plot  the  entire 
grid.  I0PTI0N  controls  the  type  of  plot.  I0PTI0N  =  0  means  a  plot 
showing  all  grids  within  the  bounding  indices.  I0PTI0N  =  1  means 
plot  velocity  vectors.  I0PTI0N  =  2  means  place  values  of  BT2  (plastic 
counter)  in  zones  which  do  not  use  equations  of  state  2  or  3  or  6. 
I0PTI0N  *  3  means  place  values  of  BT1  (PMIN  counter)  in  zones  which 
do  not  use  states  2  or  3.  XHORIZ  controls  which  variable,  X  or  Z,  will 
be  the  horizontal  axis  on  the  plots.  XHORIZ  =  0  means  the  X  axis  is 
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the  horizontal  axis.  XHORIZ  =  1  means  Z  Is  horizontal.  X6RDMIN, 
XGRDMAX,  ZGRDMIN,  ZGRDMAX  define  the  min  and  max  coordinates  for  the 
plot.  OPTION(l),  0PTI0N(2)  and  0PTI0N(3)  control  velocity  vector  plots 
They  are  not  Input  variables  In  that  they  are  automatically  set  to  0. 
They  can  be  changed  only  through  UPDATE.  Zero  values  for  these 
variables  provide  reasonable  defaults.  0PTI0N(3)  is  the  minimum 
value  of  velocity  to  display.  The  default  value  is  1.  OPTION(l) 
and  OPTION (2)  control  the  size  of  the  maximum  velocity  vector  on 
the  plot  (max  for  all  the  qrid  not  just  that  on  the  plot).  OPTION(l) 
controls  the  X  velocity  component.  0PTI0N(2)  controls  Z.  They  are 
mapping  variables  and  are  set  to  the  velocity  values  desired  for  76 
raster  points  (0.5  inches  on  the  SC4020  plotter).  Normally  the  two 
are  equal.  Independent  controls  are  provided  so  that  velocity 
vectors  can  have  their  proper  angles  even  if  XGRDMAX  -  XGRDMIN  is  not 
equal  to  ZGRDMAX  -  ZGRDMIN.  Default  values  are  such  that  the  max 
velocity  is  0.5  inches  long. 

On  cycle  0,  or  the  first  cycle  after  a  restart,  HULL  boundary 
stations  (if  any)  are  shown  on  all  grid  plots  using  the  letter  H  at 
the  station  location.  At  the  same  time  T00DY  time  history  locations 
are  shown  with  their  station  numbers. 

Figures  I I -6  and  I I -7  are  examples  of  qrid  and  velocity  vector 
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GRAPH  PLOTTING 


If  IGRAF  is  greater  than  0,  GRAPH  plots  occur  at  problem  times 
TGRFMIN,  TGRFMIN  +  TGRFDEL ,  TGRFMIN  +  2  TGRFDEL ,  ...»  TGRFMAX.  The 
user  inputs  TGRFMIN  and  TGRFDEL.  TGRFMAX  Is  set  at  1.E6.  To  change 
TGRFMAX  the  user  must  UPDATE  in  a  new  value  In  SUBROUTINE  RDGRAF. 
Each  plot,  from  1  to  IGRAF,  must  have  two  Input  cards  specified.  Th 
first  is 

IGRFMIN,  IGRFMAX,  IGRFDEL,  OGRFMIN,  JGRFMAX 
the  second  is 

OGMIN ,  OGMAX,  IORDGF 

Each  plot  will  consist  of  the  variable  specified  by  IORDGF  being 
plotted  vs  X  values  between  JGRFMIN  and  JGRFMAX  for  II  columns, 

IGRFMIN,  IGRFMIN  +  IGRFDEL,  IGRFMIN  +2IGRFDEL . IGRFMX. 

The  II  value  Is  indicated  on  each  curve.  The  ordinate  min  and  max 
values  are  OGMIN  and  OGMAX.  If  these  are  zero,  the  code  finds  min 
and  max  values  automatically.  IORDGF  values  and  the  variables 
plotted  are 

IORDGF  VARIABLE 

4  X 

5  2 

6  UX 

7  UZ 

8  TXX 

9  TYY 

10  TZZ 

11  TXZ 

12  Q 

13  EXX  (really  DQXX) 

14  EZZ  (really  DQZZ) 

15  EXZ  (really  DQXZ) 

16  RHO 

17  C  (longitudinal  sound  speed) 

18  t 

19  AL  (zone  area) 


20 

XM  (zone  mass) 

21 

BT 

22 

BT1 

23 

BT2 

24 

BT3  (normally  hoop  strain) 

25 

Total  velocity 

26 

Pressure 

Figures  I I -8  and  I I -9  are  examples  of  GRAPH  plots.  In  these 
examples,  the  routines  were  modified  to  plot  total  velocity  in  FPS 
vs  initial  radial  position,  X,and  velocity  angles  with  respect  to 
the  Z  axis  (THETA)  vs  initial  radial  position.  To  do  this  the 
array  DQXX  was  used  to  contain  initial  X  and  DQXZ  was  used  to  contain 
velocity  angle.  Then  we  simply  plotted  velocity  vs  DQXX  and  DQXZ  vs 
DQXX.  GRAPH  was  UPDATED  slightly  to  change  velocity  and  position 
units  to  FPS  and  inches.  The  X  calls  in  GRAPH  were  changed  to  DQXX 
calls.  XLAB  and  YLAB,  abscissa  and  ordinate  labels  were  changed  in 
GRAPH.  DQXX  was  set  in  SETUP  1  after  statement  395.  DQXZ  was  set 
after  statement  970  in  MAIN  by  usinq  the  call 

CALL  SUBANG  (UX(M,  INN),  UZ(M,  INN),  DQXZ(M,  INN)) 

SUBANG  is  a  built-in  utility  routine  which  computes  angle  in  degrees 
between  the  velocity  vector  and  the  positive  Z  axis.  The  angle 
returned  varies  from  -90  through  270  degrees.  Of  course,  the  coding 
in  MAIN  setting  DQXX  to  EXX  and  DQXZ  to  EXZ  had  to  be  deleted.  In 
addition  coding  after  statement  500  in  SETUP  1  had  to  be  changed  so 
that  DQXX  would  not  be  made  zero  (after  it  had  been  set  to  initial 
X  values  after  statement  395).  A  little  practice  will  allow  the  user 
to  plot  virtually  any  variable  vs  any  other  variable  which  changes 
with  JJ  values.  The  utility  Subroutine  PRINC  may  be  helpful  in  this 
regard.  It  returns  maximum  and  minimum  principal  stresses  or  strains 
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FIGURE  11-9 


given  the  XX,  ZZ  and  XZ  values.  It  can  be  called  as  follows. 

CALL  PRINC  (XX,  ZZ,  XZ,  El,  E2,  ES) 
where  XX,  ZZ  and  XZ  are  given  stress  or  strain  component  values. 

El  Is  the  maximum  stress  or  strain,  E2  Is  the  minimum  and  ES  Is 
the  maximum  shear. 

TIME  HISTORIES-AIR  FORCE 

Routines  have  been  added  to  allow  the  user  to  save  time  history 
Information  for  TOODY  zones  on  TAPE  12.  This  Information  can  then 
be  plotted  using  the  HULL  station  plotter  In  program  PULL  or  other 
plotting  routines  available  from  AFATL/DLOW.  Time  history  data  Is 
saved  In  the  same  compact  format  used  by  the  HULL  code.  This  technique 
has  been  shown  to  be  very  efficient  In  time,  memory  and  tape  require¬ 
ments  . 

Time  histories  are  controlled  by  the  Input  variable  ISPEC. 

ISPEC  should  be  set  to  the  number  of  stations  desired.  One  card 
must  be  Input  for  each  station  from  1  to  ISPEC.  The  card  contains 

IITIME,  JJTIME,  XTIME,  ZTIME 

where  (IITIME,  JJTIME)  are  the  (II,  JJ)  values  for  the  station  and 
(XTIME,  ZTIME)  are  the  stations'  (X.  Z)  values.  The  code  expects 
that  either  (II,  JJ)  or  (X,  Z)  will  be  specified.  The  unspecified 
set  Is  filled  In  by  the  code  In  SETUP!. 


The  following  data  is  saved  at  each  station  (in  the  order  shown) 


1  —  Density  with  station  number  in  the  lower  bits 

2  --  Radial  displacement 

3  --  Radial  velocity 

4  --  Axial  velocity 

5  --  Radial  acceleration 

6  --  Axial  acceleration 

7  --  Radial  stress 

8  Shear  stress 

9  --  Axial  stress 

10  --  Hoop  stress 

11  --  Radial  strain 

12  Shear  strain 

13  --  Axial  strain 

14  --  Hoop  strain 

15  --  Axial  displacement 

If  plane  geometry  is  used  (IALPH  =  1)  then  radial  components  are 
X  components.  Strains  are  computed  in  the  time  history  routines  and 
are  independent  of  zone  variables  DQXX,  DQXZ,  DQZZ  and  8T3.  These 
zone  variables  can  be  used  for  other  purposes  and  strains  can  still 
be  obtained  using  time  histories.  Displacements  are  computed  by 
assigning  initial  X  and  Z  values,  XXTO  and  ZZTO ,  to  each  station. 

The  code  assigns  these  values  as  the  node  point  values  for  (IITIME, 
JJTIME).  Displacements  are  *hen  comnuted  by  determining  nodal  X. 

7  changes.  Thus  displacements  are  ind°pendent  of  any  (XTIME,  ZTIME) 


settings.  The  (XTIME,  ZTIME)  values  are  used  only  as  position 
Indicators  on  the  plots.  If  defaulted  to  zero,  these  (XTIME,  ZTIME) 
values  are  set  to  mid-point  in  the  (IITIME,  JJTIME)  zone. 


Data  is  saved  for  any  station  whenever  total  acceleration  has 
changed  from  the  last  value  saved  by  0.2  percent.  XX  stress  is 
used  instead  of  acceleration  for  an  applied  pressure  zone.  Data  is 
packed  in  records  of  1024  words  length. 

Dimensioning  in  common  block  THIS!  assumes  no  more  than  100 
stations. 

Figures  I I -10  through  11-16  show  the  results  of  using  the  PULL 
plotter  with  a  T00DY  station  tape.  (PULL  was  modified  by  Ms  Cydney 
Westmoreland  to  plot  the  15  variables  in  the  format  shown  in  the 


figures) . 


UTILITY  R0UTINES-AIR  FORCE 


Utility  routines  SUBANG  and  PRINC  have  already  been  discussed. 

Utility  routine  FPDUMP,  written  primarily  by  Mr  Harold  Iuzzolino 
while  at  the  Air  Force  Weapons  Laboratory,  provides  decimal  dumps 
whenever  errors  occur.  It  relies  on  system  subroutines  REC0VR  and 
PDUMP.  It  is  called  whenever  a  mode  error  occurs.  Subroutine 
EREXIT  has  a  division  by  zero  built  into  it,  so  FPDUMP  is  called  also 
whenever  EREXIT  is  called.  FPDUMP  prints  out  the  location  of  the 
error,  the  values  of  all  registers  and  the  first  100  code  words, 

II,  JO,  NCYCLE  and  T  when  called  and  the  values  of  major  COMMON 

blocks  and  the  IC  columns.  As  now  set  up  it  skips  dumping  general  coding 
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between  COMMON  blocks  .To  dump  the  entire  code  change  ISHORT  to  0. 

Utility  routines  have  been  written  to  allow  zone  dropping  as 
specified  by  the  user.  These  routines  are  often  essential  to  con¬ 
tinued  running  In  very  violent  problems.  The  alternative  Is  to  use 
rezoning.  Fortunately  many  problems  become  highly  distorted  In  areas 
of  little  Interest  and  at  a  time  such  that  they  can  be  dropped  out 
of  the  grid  without  adverse  effects.  Of  course,  such  drops  have  to 
be  carefully  studied  or  they  can  actually  change  the  problem  away 
from  that  desired. 

Input  to  use  zone  dropping  routines  has  to  be  UPDATED  into 
SUBROUTINE  DROPIN  (use  *1  DR0PIT.19  followed  by  the  appropriate 
information).  The  input  information  needed  for  each  drop  is  TTDROP(K), 
I I DROP ( K) ,  JJDROP(K)  and  JTIDROP(K)  where  K  can  vary  from  1  to  20 
(that  is,  up  to  20  drops  can  occur).  TTDROP  is  the  time  of  the  drop 
in  seconds.  If  it  is  desired  to  drop  II  columns  then  JJDROP  and 

JTIDROP  must  be  0  and  IIDROP  set  to  the  number  of  II  columns  to  drop. 

Dropping  begins  with  column  1  and  goes  through  IIDROP.  Each  column 
is  entirely  dropped  from  the  problem  and  arrays  such  as  IGRDMIN  and 
IGRDMAX,  IITIME,  etc.  are  adjusted.  Coding  assumes  that  IACT  is  up 
to  IT-1  before  any  drops  occur.  Material  indicators  in  the  new  II  =  1 

column  are  changed  to  21  (free)  If  the  material  In  zone  II+l  at  the 

same  00  value  Is  real  (less  than  20).  The  indicator  is  changed  to 
that  for  II+l  at  the  same  00  value  If  the  material  indicator  at 
(II+l,  00)  Is  a  boundary  Indicator.  It  is  assumed  that  no  more  than 
one  section  of  II  lines  will  be  deleted  on  any  one  drop. 


If  it  is  desired  to  drop  JJ  rows,  then  IIDROP  is  set  to  zero 

and  JJDROP  is  set  to  the  last  real  JJ  value  it  is  desired  to  retain. 

The  drop  can  occur  in  only  one  section.  That  section  is  specified  by 
setting  JTIDROP.  For  example,  if  it  is  desired  to  drop  JJ  =  26 
through  30  (where  30  is  JMAX  in  the  JTIDROP  section)  then  set  JJDROP 
to  25.  JMAX  of  the  new  section  so  created  will  be  26.  Real  material 
will  therefore  end  at  JJ  =25.  It  is  always  assumed  that  the  drop 
goes  through  JMAX  (JTIDROP).  Arrays  such  as  JGRDMIN,  JGRDMAX,  JJTIME, 
etc.  are  automatically  updated  to  reflect  the  change.  The  boundary 
indicator  at  the  new  JMAX  is  set  to  that  at  the  old  JMAX. 

If  IACT  =  IT-1  prior  to  any  drops  then  there  will  be  no  conflicts 

if  the  problem  has  to  be  restarted.  However,  the  restart  deck  should 
have  new  values  for  IT,  JT,  GRID  and  GRAPH  indices  and  section 
definition  (the  card  beginning  with  NJVAR).  Material  indicator  cards 
and  geometry  definition  cards  need  not  be  changed  from  their  old  values 
The  history  station  indices  can  also  remain  at  their  old  values. 

Subroutine  DFPRINT  places  a  message  in  the  dayfile  whenever  a 
restart  file  is  made,  the  problem  terminates  or  SWITCH  2  is  ON.  The  ■ 
message  gives  IACT.NCYCLE,  T  and  DT  values  plus  NFILE  (the  restart 
file  being  made) . 

PPCALL  is  a  local  ADTC  system  routine.  As  used  in  TOODY  it 
returns  the  job  card  time  limit  so  that  TOODY  can  exit  gracefully 
(i.e.  by  ending  files,  doing  plotting  chores,  etc.)  prior  to  CP  time 


CODE  TERMINATION-AIR  FORCE  AND  SANDIA 


TOODY  will  terminate  If  T  Is  greater  than  or  equal  to  TMAX, 
SWITCH  1  Is  ON,  job  card  time  limit  is  only  30  seconds  away  or  If 
CP  time  exceeds  RSMAX.  At  termination,  a  restart  file  Is  made  If 
r.SMIN  is  non-negative,  GRID  and  GRAPH  plotting  are  accomplished  If 
called  for  by  Input,  printing  Is  accomplished  If  called  for  by  Input 
and  time  station  buffer  areas  are  dumped  on  TAPE  12. 

INTERMEDIATE  VARIABLE  STORAGE-AIR  FORCE 


Intermediate  variables  are  stored  only  on  disc  (TAPE  21  and 
TAPE  22).  Some  versions  of  TOODY  use  Extended  Core  Storage  (ECS). 
Since  the  AFATL  has  no  ECS,  this  version  of  TOODY  relies  completely 
on  disc.  It  would  be  fairly  simple  to  modify  the  code  to  use  ECS. 
BUFFER  IN  and  OUT  statements  could  be  changed  to  READEC  and  WRITEC 
statements.  All  IF  UNIT  checks  on  NTP1  and  NTP2  would  have  to  be 
deleted. 

As  now  written,  the  AFATL  TOODY  uses  overlapping  buffering  and 
numerical  computing  -  l.e.,  a  buffer  operation  Is  started  and  the 
code  does  not  wait  for  It  to  be  finished  (an  IF  UNIT  check)  until 
the  data  Is  absolutely  required  for  further  computing. 

BOWTIE  DISTORTION  CORRECTION -SANDIA 


TOODY  has  a  built  In  viscosity  term  to  dampen  hourglass  or 
bowtle  zone  distortions.  The  viscosity  is  controlled  by  IKEYSW 
(1  for  on,  0  for  off)  and  XKCON  (the  amount  of  viscosity).  IKEYSW 
Is  set  to  1  and  XKCON  to  0.25  automatically.  XKCON  (set  in  SETUP) 
should  never  exceed  1.  (XKCON  Is  multiplied  by  0.1  in  MAIN  prior  to 
use) . 
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DEBONDING/REBONDING  SLIDE  LINE-SANDIA 

A  slide  line  option  has  been  added  by  Sandia  to  allow  slide  lines 
to  collide,  debond  and  rebond.  The  material  indicator  for  this  type 
of  line  is  MAT=26.  Collision  and  separation  must  proceed  only  from 
points  bounding  the  area  In  contact.  Separation  proceeds  one  left 
hand  mesh  per  cycle.  Collision  proceeds  until  the  first  non-collision 
is  found.  Additional  input  is  required  to  use  this  option.  The  arrays 
JMINL,  JMAXL ,  JMINR  and  JMAXR  must  be  input  for  each  interface.  In 
addition,  separation  stresses  must  be  input  in  the  array  SIGSEP. 

JMINL  and  JMAXL  define  the  JJ  values  away  from  the  region  of  contact 
for  the  left  side  of  the  slide  line.  JMINR  and  JKAXR  provide  the  same 
definition  for  the  right  side.  Debonding  occurs  when  normal  stress 
across  the  slide  line  exceeds  the  SIGSEP  value  for  that  line. 

MOMENTUM  CALCULATION-SANDIA 

Sandia  has  added  a  more  accurate  momentum  zone  calculation 
than  specified  in  References  1  and  2.  It  no  longer  assumes  that  the 
areas  and  masses  of  the  momentum  zone  are  one  fourth  of  the  totals 
for  the  four  zones.  Instead  it  computes  the  area  and  mass  each 
time  they  are  required.  This  improves  the  accuracy  of  the  momentum 
calculation  for  non-rectangular  zoning  or  rectangular  zoning  which  has 
been  distorted. 

JACT-AIR  FORCE 

John  Levesque,  while  at  AFWL ,  added  coding  for  an  activity  counter 


along  JJ  rows.  The  activity  counter  is  controlled  by  the  input 
variable  JACT.  Each  II  column  has  its  own  activity  counter,  J KACT (II). 


Initially,  all  JKACT  values  are  set  to  JACT.  After  that  they  each 
advance  by  comparing  pressure  at  JKACT (II )+l  to  PACT.  Computations 
are  made  along  each  II  column  through  JKACT ( 1 1 )+l .  Zone  variables 
from  JKACT (II )+2  through  JT  are  kept  at  their  old  values.  This 
option  has  not  been  thoroughly  checked  in  AFATL  TOODY.  There  may  be 
some  problem  in  using  it  with  other  features.  Presently,  JACT 
is  always  set  to  0,  bypassing  the  option.  If  the  user  desires  to 
calculate  with  JACT  he  can  UPDATE  a  non-zero  value  into  SETUP. 

PRINT  EDITS-AIR  FORCE 

Printing  is  controlled  by  IPRINT  as  in  References  1  and  2. 

In  addition,  the  code  now  prints  out  the  weight  of  each  material 
in  grams  and  pounds  at  SETUP  time.  TOODY  uses  one  radian  of  arc 
around  an  axis  of  rotation.  The  weight  printout,  however,  multiplies 
by  2H  so  that  the  complete  weight  of  the  object  being  modelled  is 
given.  It  also  prints  out  mass  averaged  velocity  components  every 
time  GRIDS  is  called. 

UNITS-AIR  FORCE 

It  is  highly  recommended  that  only  CGS  units  be  used  as  inputs. 
These  are: 

Dimensions  cm 

Stress  dynes/cm 


Energy  Density  ergs/gm 


There  are  some  routines  in  T000Y,  such  as  the  initial  weight 
computations,  which  assume  CGS  is  used. 

INPUT  CARD  FORMAT- AIR  FORCE 

Input  is  now  format  free  and  uses  the  READ  *  option  in  FTN. 

Under  this  option  input  can  be  in  any  fields  on  the  card.  Input 
variables  are  separated  by  blanks  or  commas.  Floating  point 
variables  can  be  input  as  integers,  numbers  with  a  decimal  point 
or  in  E  format.  For  example,  the  number  100  can  be  input  as  100,  100. 
or  1.E2.  Integers  must  not  be  E  type  or  have  decimal  points  in  their 
field. 

All  numbers  to  be  read  by  a  READ  *  command  must  be  on  the 
card.  Blank  fields  are  not  understood  since  blanks  are  delimiters. 

PREPROCESSOR* SANDIA  AND  AIR  FORCE 

Sandla  added  a  preprocessor  program  to  read  TOODY  Input  and  size 
some  rrays,  particularly  the  IC  array.  This  preprocessor  has  been 
extended  by  the  Air  Force  and  deletes  unnecessary  equations  of  state, 
sliding  routines  if  there  are  no  slip  surfaces  in  the  problem,  etc., 
in  addition  to  sizing  many  arrays.  The  preprocessor  makes  a  file 
named  DCJID  which  contains  UPDATE  commands.  TOODY  uses  this  file 
through  a  *READ  DCJID  card  in  the  TOODY  UPDATE  deck.  This  system 
insures  that  TOODY  is  no  larger  than  required  for  any  given  problem. 


SECTION  III 


INPUT 

AFATL  TOODY  input  is  discussed  below.  Unless  specified,  input 
is  format  free. 

CARD  1  (A10  Format) 

Problem  ID  line  -  blank  If  a  restart 
CARD  2 

IT,  JT,  NOMAT,  IALPH,  I GRID,  IGRAF,  I PRINT,  ISPEC,  IACT 

where 

IT  is  the  max  II  value  in  the  problem 
JT  is  the  max  JJ  value  if  sections  are  not  used.  JT  *  0 
if  sections  are  used. 

NOMAT  is  the  number  of  real  materials 
IALPH  is  the  symmetry  Indicator. 

0  if  the  X  axis  is  an  axis  of  rotation 

1  if  plane  geometry 

2  if  the  Z  axis  is  an  axis  of  rotation 
IGRID  is  the  number  of  GRID  type  plots 
IGRAF  is  the  number  of  GRAPH  plots 

IPRINT  is  1  for  print  output,  0  otherwise. 

ISPEC  is  the  number  of  time  history  stations 

IACT  Is  the  Initial  value  of  the  II  activity  counter. 

If  IACT  =  0,  it  is  set  to  IT-1.  It  is  advanced  by 
comparing  pressure  in  the  zones  along  IACT+1  with  PACT. 
PACT  is  set  automatically  to  100  dynes/cm^. 


TMAX,  I HULL  IN ,  RHOFF,  YHOFF 

inhere 

TMAX  is  max  problem  time  in  seconds. 

IHULLIN  is  1  if  HULL  input  is  to  be  used. 

RHOFF,  YHOFF  are  radial  and  axial  offsets  to  be  added  to  HULL 
station  coordinates. 

NOTE:  There  are  many  variables  which  would  be  input  at  this  time 

in  the  code  specified  by  References  1  and  2,  but  which  are  now  set 
automatically.  To  change  these  values  one  would  have  to  UPDATE  in 
new  values  in  SETUP.  These  are 

PMAX  =  50. El 2  (max  pressure  allowed) 

XKX  =  0.8  (time  step  multiplier) 

B1  *  1.2  (bulk  quadratic  viscosity  coefficient) 

B2  =  0.6  (bulk  linear  viscosity  coefficient) 

B3  =  B4  =  0  (deviator  viscosity  coefficients) 

DELTMIN  =  l.E-12  (min  time  step  allowed) 

DELT  -  1.E6  (code  then  will  choose  initial  time  step) 

DDT  =1.1  (time  step  max  growth  factor) 

IECSW  =  0  (enerq.y  check  switch) 

IRRSSW  =  0  if  ID  line  not  blank 
=  1  if  ID  line  blank 
IRSF  =  1  (restart  file  number) 
iuNORDT  =  0 

JACT  =  0  (this  bypasses  all  JJ  activity  checks) 
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CARD  3  (Present  only  if  I PRINT  =  1 


TPDEL ,  IPMIN,  IPMAX,  IPDEL,  JPMIN,  JPMAX,  JPDEL 

where 

TPDEL  is  the  increment  In  problem  time  at  which  print  edits 
are  desired. 

Edits  are  desired  for  II  values,  IPMIN,  IPMIN  +  IPDEL,  . ..,  IPMAX 

and  JJ  values  JPMIN,  JPMIN  +  JPDEL,  ....  JPMAX 

CARD  4  (present  only  if  ISPEC  greater  than  0 

NOTE:  There  are  ISPEC  numbers  of  Card  4.  Each  card  applies  to  a 

single  station. 

IITIME,  JJTIME ,  XTIME ,  ZTIME 

where 

IITIME,  JJTIME  is  the  II,  JJ  for  the  time  history  station. 

XTIME,  ZTIME  is  the  X.Z,  for  the  station. 

Either  II,  JJ  or  X,  Z  can  be  specified  with  zero  values  for  the  other 
unspecified  set. 

CARD  5  (present  only  if  IGRID  greater  than  0) 

TGRDDEL 

where 

TGRDDEL  is  the  problem  time  Increment  at  which  GRID  plots 
are  desired. 

CARDS  6  and  7  (present  if  IGRID  greater  than  0) 

NOTE:  There  are  IGRID  sets  of  Cards  6  and  7 


IGRDMIN,  IGRDMAX,  JGRDMIN,  JGRDMAX,  IOPTION 
XGRDMIN,  XGRDMAX,  ZGRDMIN,  ZGRDMAX ,  XHORIZ 

When 

IGRDMIN  is  the  minimum  II  value  for  the  plot 
IGRDMAX  is  the  maximum  II  value  for  the  plot 
JGRDMIN  is  the  minimum  JJ  value  for  the  plot 
JGRDMAX  is  the  maximum  JJ  value  for  the  plot 
IOPTION  is  0  for  grid  plot 

1  for  velocity  vector  plot 

2  for  a  plot  showing  zones  which  have  ever  been  plastic 

3  for  a  plot  showing  zones  which  have  ever  had  a 

pressure  less  than  PMIN 

XGRDMIN,  XGRDMAX  are  min  and  max  values  of  X  desired  on  the  plot 

ZGRDMIN,  ZGRDMAX  are  min  and  max  values  of  Z  desired  on  the  plot 

XHORIZ  is  0  if  the  X  axis  should  be  horizontal  on  the  plot. 

It  is  1  if  the  Z  axis  should  be  horizontal. 

CARD  8  (present  only  if  I GRAF  greater  than  0) 

TGRFMIN,  TGRFDEL 

where 

TGRFMIN  is  the  first  problem  time  at  which  GRAPH  plots  are  desired 
TGRFDEL  is  the  problem  time  increment  for  plots. 

CARDS  9  and  1 0  (present  if  IGRAF  greater  th an  0) 

There  are  IGRAF  sets  of  cards  9  and  10  required. 


IGRFMIN,  IGRFMAX ,  IGRFDEL  JGRFMIN .  JGRFMAX 
OGMIN,  OGMAX ,  IORDGF 


where 


IGRFMIN  Is  min  II  for  the  plot 
IGRFMAX  is  max  II  for  the  plot 
IGRFDEL  is  the  increment  in  II 
JGRFMIN  is  min  JJ  for  the  plot 
JGRFMAX  is  max  JJ  for  the  plot 
OGMIN  is  min  value  for  the  ordinate 
OGMAX  is  max  value  for  the  ordinate 
If  OGMIN®OGMAX*0  then  the  code  assigns  min  and  max  values 
IORDGF  indicates  which  zone  variable  to  plot  vs  X.  The 
possible  input  values  are 


IORDGF 

VARIABLE 

4 

X 

5 

Z 

6 

UX 

7 

uz 

8 

TXX 

9 

TYY 

10 

TZZ 

11 

TXZ 

12 

Q 

13 

EXX 

(really  DQXX) 

14 

EZZ 

(really  DQZZ) 

15 

EXZ 

(really  DQXZ) 

16 

RHO 

17 

C 

18 

E 

19 

AL 

20 

XM 

21 

BT 

22 

BT1 

23 

BT2 

24 

EYY(really  BT3) 

25 

Total  Velocity 

26 

Pressure 
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CARD  11  (present  only  If  IECSM  has  been  UPDATED  to  a  value  of  1) 

EPE,  ECON,  NSCYCLE 

where 

EPE  is  the  relative  energy  error  allowed. 

ECON  is  the  value  of  energy  below  which  no  check  will  be  made. 
NSCYCLE  is  the  cycle  to  begin  checking 
NOTE:  Use  of  this  energy  check  option  is  usually  not  required.  If 
instabilities  occur,  they  usually  grow  very  rapidly  resulting  in 
violation  of  the  PMAX  criterion  within  a  few  cycles.  There  are  other 
difficulties  with  using  an  energy  check  such  as  the  fact  that  applied 
pressure  work  is  not  considered. 

CARD  12 

RSMIN,  RSMAX,  RSDEL 

where 

RSMIN  is  any  negative  number  such  as  -1  if  restart  files 
are  not  desired 

is  problem  time  for  first  dump  if  RSMAX  is  less  than  0 
is  CP  running  time  for  first  dump  if  RSMAX  greater  than  0 
RSMAX  is  running  time  in  CP  minutes.  Use  a  negative  value  if 
restart  dumps  are  desired  at  certain  problem  times. 

In  this  case  the  absolute  value  of  RSMAX  is  used  for  CP 
time  limit. 

RSDEL  is  the  increment  in  problem  time  at  which  restart  dumps 
are  desired  if  RSMAX  is  less  than  0. 

is  the  increment  in  CP  time  (minutes)  if  RSMAX  is  greater 


than  0. 
Some  examples  are 


I 


-1,  30,  0  would  mean  no  restart  dumps  desired  and  a  max 
CP  running  time  of  30  minutes.  (In  this  case  RSDEL  has  no 
meaning  and  can  be  set  to  0) 

10.E-6,  -30,  20.E-6  would  mean  dumps  are  desired  every 
20.E-6  seconds  problem  time  beginning  at  10.E-6  seconds  and  that 
30  minutes  Is  the  CP  time  limit  desired. 

10,  120,  20  would  mean  dumps  are  desired  after  every 
elapsed  20  minutes  CP  time  beginning  at  10  minutes.  The 
problem  should  run  no  more  than  120  minutes. 

CARD  13  (present  only  if  JT  =  0  on  Card  2) 

NJVAR,  IJ(1),  JMIN(l),  JMAX(1 ),...,  IJ(NJVAR),  JMIN(NJVAR), 
JMAX(NJVAR) 

where 

NJVAR  Is  the  number  of  sections  in  the  problem 
10(1)  is  the  last  II  value  in  the  first  section 
JMIN(l)  is  the  min  JO  value  in  the  first  section 
JMAX  (1)  is  the  max  00  value  in  the  first  section 


n 


IJ(NJVAR)  is  the  last  II  value  in  the  last  section  (IT)  ? 

JMIN(NJVAR)  is  the  min  JJ  value  for  the  last  section 
JMAX(NJVAR)  is  the  max  JJ  value  for  the  last  section 

» 
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NOTES:  Up  to  10  sections  are  allowed. 

It  is  generally  better  to  set  JMIN  values  to  1  if  at  all 
possible. 

CARD  14  (there  can  be  up  to  40  of  the  following  material  definition 


cards' 


KM,  IS1,  IS2,  JS1,  JS2 

where 

KM  is  material  number,  including  boundary  indicators 

151  is  the  starting  II  value  for  KM 

152  is  the  ending  II  value  for  KM 
JS1  is  the  starting  JJ  value  for  KM 
JS2  is  the  ending  JJ  value  for  KM 

CARD  15  (Material  definition  termination  card) 


This  card  terminates  read  in  of  material  definition  input.  It  must 
be  one  of  the  following.  Either 

0  0  0  0  0 


0-1  000 

The  latter  is  used  if  geometry  definition  cards  are  to  be  input  later. 
The  former  is  used  if  material  definition  cards  also  determine  geometry 
definition. 

CARD  16  (present  if  KM  =  25  or  26  or  any  material  definition  cards) 
FRICT(l),  FRICT(2) .  FRICT(ISLIDE) 


where 


FRICT(l)  is  the  Coulomb  friction  coefficient  for  slide  linel 
(i.e.,  the  slide  line  with  the  smallest  II  value) 
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FRICT(2)  Is  the  coefficient  for  the  second  slide  line. 


FRICT  (ISLIDE)  is  the  coefficient  for  the  last  slide  line  In 
the  problem. 

i'JTES :  There  can  be  no  more  than  8  slide  lines. 

FRICT  -  -1  means  infinite  friction  (no  sliding) 

■  0  means  perfect  no  friction,  sliding 
>  0  means  slide  with  the  friction  coefficient  specified 


CARD  16A  (present  if  KM  -  26  on  any  material  definition  cards) 
JMINL(l) 

JMAXL(l) 

JMINR(l) 

JMAXR(l) 


JMINL(8) 

JMAXL(R) 

JMINR (8) 

JMAXR(8) 

CARD  16B  (present  if  KM  *  26  on  any  material  definition  cards) 
SIGSEP(I) 

SI GSEP ( 2 ) 


SIGSEP(8) 
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resent  only  if  the  terminating  material  Indicator  card  was 


0.  zh  o.  o.  Q) 

KG,  ISG1 ,  ISG2 ,  JSG1 ,  JSG2 

where 

KG  is  geometry  region  number  (beginning  with  1). 

ISG1  is  min  II  value  for  region  KG 

ISG2  is  max  II  value  for  region  KG 

JSG1  is  min  JJ  value  for  region  KG 

JSG2  is  max  JJ  value  for  region  KG 

CARD  18  (Geometry  termination  card.  Present  only  if  CARDS  17  present) 
This  card  must  be  five  zeros. 

0  0  0  0  0 

The  next  set  of  cards  define  equation  of  state  properties  for  material 
1  through  N0MAT.  There  are  varying  numbers  of  cards  depending  on  the 
equation  of  state  option.  There  must  be  N0MAT  sets  of  these  cards. 
CARD  19 
IES 

where 

IES  is  equation  of  state  indicator 

1  for  state  1  (elastic/plastic) 

2  for  state  2  (high  explosive) 

3  for  state  3  (ideal  gas) 

4  for  state  4  (earth  media) 

(5  is  not  used) 

6  for  state  6  (metal  foam) 

7  for  state  7  (GRAY  metal) 

9  for  state  9  (slurry  EOS) 


S  Is  the  slope  of  the  US/UP  Hugonlot,  If  S  =  0, 
PH  calculated  as  poco2y 
ro  is  initial  Gruneisen  ratio 


NOY  Is  strength  indicator 
0  -  no  strength 

1  -  elastic/perfectly  plastic 

2  -  perfectly  elastic 

3  -  energy  variable  plastic 

4  -  linear  work  hardening 

5  -  nonlinear  work  hardening 


YO  is  initial  plastic  flow  stress 


where 

Y1  is  Young's  modulus 
Y2  is  the  tangent  modulus 

NOTES:  This  input  has  been  greatly  simplified  from  Reference  1  by 

making  certain  assumptions.  These  are 

1.  Rigidity  modulus,  6,  is  always 


where  K  is  current  bulk  modulus. 

* 

2.  Gruneisen  ratio,  F  is  always 

r  =  ropo/p 

3.  PH  is  either 

PH  =  pocp^n/(l  -  Sn)^ 
or 

2 

PH  =  poco  ij 

4.  TDEP  always  -1  (no  energy  deposition) 

5.  EH  always  1.E1Q0  (no  melt  energy  changes  on  PMIN  or  Y) 

6.  ESUB  always  l.Einn  (no  switch  to  a  vapor  equation  of  state  at 
high  energys) 

If  these  built-in  defaults  are  not  desired,  the  READ  cards  can  be 
changed  in  SETUP  or  changes  in  equation  of  state  values  can  be 
UPDATED  into  the  code.  If  changes  are  desired,  they  should  be  made  in 
the  ESC  array  for  the  material.  This  array  is  as  follows  for  material  K 


ESC  (  K,  1)  =  po 
ESC  (K,  2)  =  co 

ESC  (K,  3)  =  QDEP  (energy/unit  mass  to  deposit) 
ESC  (K,  4)  =  PMIN 

ESC  (K,  5)  =  TDEP  (energy  deposition  time) 

ESC  (K,  6)  =v 

ESC  (K,  7)  not  used 

ESC  (K ,  8)  *  NOK 

If  NOK  =  0  then  PH  =  poco2n/0  -  Sn)2 
if  NOK  between  1  and  6,  then 
PH  =  poco^O  +  k-jn  +  K2n2  +  *•*  +  KN0KnN0K) 
if  NOK  =  7,  then 

PH  =  poco2y(l  +  K-|U  +  K2u2  +  *••) 

ESC  (K,  9)  *  poco2  (filled  In  by  the  code) 

ESC  (K,  10)  =  S  if  NOK  =  0 
=  If  NOK  >  1 
ESC  (K,  11)  =  K2 
« 

ESC  (K,  14)  =  K4 
ESC  (K,  15)  *  NOH 

If  NOH  =  0,  r  =  To 
If  NOH  >  1 , 

r  =  rod  +  H1  n  +  H2n2  +  •••  +  HN0HnN0H) 
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ESC  (K,  32)  =  Y2 

ESC  (K,  33)  =  Y3  (used  only  for  NOY  =  5) 

ESC  (K,  34)  =  EMELT 

ESC  (K,  35)  =  ESUB 

ESC  (K,  36)  =  y-1  for  vapor  equation 
"  EMELT,  energy  density  at  melt,  is  specified  then  Y  and  PMIN  are 
:iul  tipi  led  by 

1  -  E/ EMELT 

prior  to  use. 

ESUB  is  the  vaporization  energy  density  and  is  used  in  the  vapor 
equation  with  y-1  (See  Reference  1) 

If  IES  =  2,  then  Input. 

CARD  19A 

po,  D,  D,  y,  0.  0,  0,  0  if  gamma-law  explosive  desired. 

D  is  detonation  velocity.  Energy  of  detonation  is  computed 
internally  as  D2/ (2(y^  -  1) 
po,  D,  E,  A,  B,  R1 ,  R2,  W  If  JWL  explosive  desired. 

CARD  19B 

ND,  XD,  ZD,  TOBURN,  K0 
ND  Is  detonation  option 

ND  a  1  ignite  from  line  Z  =  ZD 
ND  a  2  ignite  from  line  X  =  XD 
ND  a  3  Ignite  from  point  (XD,  ZD) 

ND  *  4  ignite  all  explosive  at  the  same  time. 
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TOBURN  Time  for  explosive  to  begin  igniting  (normally  0) 

KO  Bulk  modulus  for  unignited  explosive.  If  non-zero 

pressure  computed  from  P  =  Koy 
If  Kc  =  0,  pressure  computed  from 
P  =  poco^n/O  -  sn 

where  po,  co  and  S  are  appropriate  for  composition  B. 
If  IES  =  3,  then  Input 

CARD  19A 

po,  co ,  y,  0,  0,  0,  0,  0 

If  IES  =  4,  then  Input 

CARD  19A 

po,  Cp,  PMIN , v  ,  PELAS,  PLOCK,  pLOCK 

where 

po  =  initial  density 
Cp  =  initial  longitudinal  wave  velocity 
PMIN  =  minimum  pressure  allowed 
PELAS  =  max  pressure  prior  to  crushing 

PL0CK,  pLOCK  =  define  intersection  of  crush  curve  with  lockup 
curve 

CARD  19B 

CoL,  S,  Y0,  YV,  PV,  GMAX,  PMID,  pM I D 

where 


CQL  and  S  define  the  lock-up  Hugoniot,  U$  =  C0|_  +  SUp 


YO,  YV,  PV  define  the  yield  surface 

GMAX  is  limiting  rigidity  modulus  (input  0  if  no  limit  desired) 
PMID,  yMID  =  defines  crushup  curve  between  PELAS,  pELAS  and  PLOCK, 
uLOCK 

NOTES: 

To  run  as  a  fluid  medium  input  v  =  0.5,  YO  =  YV  s  PV  =  0 
2.  To  u<e  the  built-in  concrete  equation  of  state,  the  first  card, 

CARD  19A,  should  be  all  O' s.  In  this  case,  CARD  19B  is  not  required. 

If  IES  =  6,  then  Input 

CARD  19A 

pos  C0,  QDEP,  PMIN ,  TDEP,  B5,  B6 

where 

pos  =  initial  density  of  the  solid 
C0  *  initial  bulk  sound  speed  of  solid 
QDEP  =  energy  density  to  be  deposited 
PMIN  =  minimum  pressure 

TDEP  =  energy  deposition  time  (-1  if  no  deposition) 

B5,  B6  are  extra  viscosity  terms. 

B5  =  3,  B6  =  0.7  are  good  values 
M D  19B 


NOK,  0,  K  or  S,  K2,  K3,  K4,  K5 

Hugoniot  definition  for  the  solid 


■ 


« 


CARD  19C 

NOH,  To,  HI,  H2,  H3,  H4,  H5 
Gruneisen  ratio  definition 
CARD  19D 

Ce»  »  ^eo’  ^so’  ^ >  E0,  0 

where 

Ce  =  initial  foam  sound  speed 
oq  =  initial  distention  ratio 

Peo  =  elastic  limit  pressure  for  unheated  material 

Ps0  =  pressure  at  which  all  voids  are  filled  for  unheated  material 

K  =  coefficient  in  energy  dependence  equations  for  Pso  and  Pe0. 

Set  to  -2  if  it  is  input  as  a  0. 

E0  =  energy  at  which  foam  strength  becomes  negligible.  Set  to 
0  for  energy  independent  model. 

CARD  19E 

ESUB,  y-1,  0,  0,  0,  0,  0,  0 
Terms  for  vapor  equation 

If  IES  =  7,  then  Input 

CARD  19A 

po,  Co,  QDEP,  PMIN,  TDEP ,  ,  0 

where  TDEP  =  -1  for  no  energy  deposition 
CARD  19B 

S,  yo,  a,  GE,  TMO,  AW,  RJ 


a 

f 

ft 


» 

¥ 


V 


I 


t 


CARD  19C 


RB,  AY,  TH,  TRaAG,  0,  0,  0 
CARD  19D 

NOG,  0,  gr  g2,  g3,  g4,  g5 
CARD  19E 

NOY,  YO,  Yl,  Y2,  Y3,  0,  0,  0 
The  ESC  array  for  State  7  is  as  follows: 


ESC 

(K, 

1)  = 

po 

ESC 

(K, 

2)  = 

Co 

ESC 

(K, 

3)  = 

QDEP 

ESC 

(K, 

4)  = 

PMIN 

ESC 

(K, 

5)  = 

TDEP 

ESC 

(K, 

6)  = 

V 

ESC 

(K, 

7)  = 

Not  Used 

ESC 

(K, 

8)  = 

S 

ESC 

(K, 

9)  = 

Yo 

ESC 

(K, 

10) 

=  a 

ESC 

(K, 

ID 

*  GE 

ESC 

(K, 

12) 

=  TMO 

ESC 

(K, 

13) 

=  AW 

ESC 

(K, 

14) 

=  RJ 

ESC 

(K, 

15) 

=  RB 

ESC 

(K, 

16) 

=  AY 

ESC 

(K, 

17) 

=  TH 
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ESC  (K,  18)  =  TRFLAG 

ESC  (K,  19)  =  Not  Used 

ESC  (K,  20)  =  Not  Used 

ESC  (K,  21)  =  Not  Used 

ESC  (K,  22)  =  NOG 


ESC  (K,  23)  =  GO  (Code  fills  in) 
ESC  (K,  24)  =  g., 

ESC  (K,  25)  =  g2 

ESC  (K,  26)  =  g3 

ESC  (K,  27)  =  q4 

ESC  (K,  28)  -  g5 

ESC  (K,  29)  =  NOY 

ESC  (K,  30)  =  YO 

ESC  (K,  31)  *  Y1 

ESC  (K,  32)  =  Y2 

ESC  (K,  33)  =  Y3 

ESC  (K ,  34)  =  Not  Used 

ESC  (K,  35)  «  Not  Used 

ESC  (K,  36)  =  Not  Used 


CARD  19A 


If  IES  =  9,  then  Input 

po,  D,  E,  AT,  Bl,  Rll ,  R21 ,  W1 


where  D  is  detonation  velocity,  E  is  energy  per  unit  mass  released 
upon  initial  burn,  A1  through  Vi!  are  parameters  for  the  first  JWL  EOS 


:;ARD  19B 


TRAMP,  EZDEL ,  A2,  B2,  R21 ,  R22s  W2 


Where  TRAMP  is  the  time  required  to  release  the  delayed  energy 
EZDEL  and  A2  through  W2  are  parameters  defining  the  second  JWL 
equation. 

CARD  19C 

ND,  XD,  ZD,  TOBURN,  B5 


where  ND,  XD,  ZD  and  TOBURN  have  the  same  meaning  as  in  State  2. 

B5  is  a  control  on  time  required  to  completely  burn  a  zone.  It 
should  be  input  as  2.5. 

CARDS  20,  20A  and  20B  are  required  to  define  shape  definition  for 
each  geometry  region.  The  cards  are  arranged  in  order  for  geometry 
regions  1  through  N6E0M  (the  last  geometry  region) 

CARD  20 
I  SET 


where  ISET  =  shape  option.  This  is  an  Integer  from  1  through  10. 
In  addition,  the  following  two  cards  are  required. 

CARD  20A 


X0  XC0N(1)  XC0N( 2 )  XC0N(3)  XC0N(4)  JMIN  JMAX 


ZO  ZCON(l)  ZCON (2)  ZCON(3)  ZC0N(4)  IMIN  IMAX 
where  IMIN,  JMIN,  IMAX  and  JMAX  define  the  boundaries  of  the  region 
NOTE:  Velocities  are  no  longer  input  but  must  be  added  using 

UPDATE  if  non-zero  values  are  required. 

If  ISET  =  7  a  third  card  is  required  as  shown  below. 

CARD  20C 

R12,  R23,  R34,  R41 ,  e^j ,  e^j 

The  specific  meaning  of  XO,  XC0N(1),  etc.,  varies  with  shape  option 
Section  II  of  this  report  and  Reference  1  explain  these  meanings. 


SECTION  IV 


CONTROL  CARDS 

The  control  card  setup  for  running  under  SCOPE  3.4  Is  as  follows: 
JOB  CARD 

CARDS  NECESSARY  TO  BRING  IN  TDPREP  and  TDMAIN  OLDPL  files. 
COPYCR,  INPUT,  TAPE  50. 

REWIND,  TAPE  50. 

UPDATE,  F,  P  =  TDPREP,  W. 

FIN,  A,  I  =  COMPILE. 

LGO. 

REWIND,  LGO. 

UPDATE,  F,  P  =  TDMAIN,  W. 

FIN,  A,  I  =  COMPILE 

CARDS  NECESSARY  TO  ASSIGN  ANY  REAL  TAPES  TO  TAPE  12, 

TAPE  9,  TAPE  25,  TAPE  26  and/or  TAPE  27 

LGO. 

7/8/9 

TOODY  INPUT  DECK 
7/8/9 

*  COMPILE  PREPRS 

ANY  PREPROCESSOR  UPDATE  CARDS 
7/8/9 

*  COMPILE  TOODY  3 

*  READ  DCJID 

ANY  TOODY  MAIN  UPDATE  CARDS 
6/7/8/9 


SECTION  V 


SOME  TIPS  ON  SUCCESSFUL  RUNNING 

A  number  of  lessions  have  been  learned  by  many  users  over  the 
years.  This  section  discusses  some  of  the  more  important  of  these. 

At  least  two  zones  are  required  across  any  material  region. 

A  single  zone  can  cause  numerical  difficulties  and  return  an  Incorrect 
answer.  Two  zones  is  the  absolute  minimum  and  can  be  used  only  in 
cases  in  which  the  material's  mass  is  the  dominant  feature--as  for 
example  in  bomb  casing  being  driven  off  by  explosives.  Three  to 
five  zones  are  needed  for  any  reasonable  shock  definition  across 
the  material. 

Zone  aspect  ratios  {length  to  width)  should  not  exceed  5. 

The  peak  pressure  in  an  explosive  zone  will  be  a  function  of 
its  size.  Zones  of  size  0.1cm  x  0.1cm  are  required  to  achieve  90% 
or  so  of  the  Chapman  -  Jouget  pressure.  On  the  other  hand,  peak 
pressure  is  often  less  important  than  impulse.  In  this  case  zones 
of  0.2  to  0.25cm  widths  and  lengths  are  adequate. 

The  JWL  explosive  formulation  should  be  used  whenever  possible. 

A  constant  gamma  explosive  will  impart  20  to  30  percent  too  much 
impulse. 

If  there  is  two  large  a  jump  in  the  factor  me  (mass  times  sound 
speed)  across  a  material  interface,  the  wrong  answer  can  result. 


Ideally,  me  changes  should  be  limited  to  no  more  than  a  factor  of  5. 

A  frictionless  slide  surface  between  materials  decouples  the  momentum 
eqii  on  on  both  sides  so  that  this  me  rule  can  be  violated  at  that 
location. 

The  accuracy  of  the  momentum  equation  approximation  deteriorates 
as  zones  become  highly  distorted.  Such  zones  should  be  dropped  if 
possible  or  rezoned  to  more  regular  shapes. 

State  1  in  TOODY  was  designed  for  stresses  up  to  500kb  or  so  and 
temperatures  up  to  one-half  of  the  material's  melt  temperature.  Changes 
in  Y  and  PMIN  with  energy  and  the  vapor  equation  of  state  are  inaccurate 
in  all  but  their  most  gross  features.  If  material  is  expected  to 
melt  or  vaporize  better  equations  should  be  used.  Such  equations  have 
been  developed  by  Sandia  (the  ANEOS  routines)  and  by  Lawrence  Livermore 
Laboratory  (the  GRAY  routines)  for  several  metals. 

It  is  recommended  that  PMIN  be  set  so  low  for  metals  that  it  is 
never  reached  --  e.g.,  at  -1.E13  or  so.  PMIN  is  not  really  a  good 
way  to  handle  spall.  If  spall  cannot  be  handled  with  a  debonding 
slide  surface,  it  cannot  be  handled  correctly.  Short  of  a  slide 
surface,  the  best  approach  is  to  monitor  principal  stress  and  strain 
in  a  zone  and  decide  if  the  zone  has  fractured  based  on  considering 
both  values.  The  momentum  routine  (starting  with  Statement  800 
in  MAIN)  can  then  be  changed  to  use  zero  stresses  when  the  zone  is 
fractured.  In  this  way,  neighboring  zones  will  encounter  a  free 


surface  (zero  stresses)  and  at  least  they  will  move  correctly. 

Such  changes  can  be  UPDATED  into  TOODY  very  easily. 

To  avoid  negative  areas  a  vector  in  the  direction  of  increasing 
II  crossed  into  a  vector  in  the  direction  of  increasing  JJ  must  have 
the  same  sign  as  a  unit  vector  in  the  Z  direction  crossed  into  a 
unit  vector  in  the  X  direction.  This  results  In  most  users  automati¬ 
cally  letting  Z  increase  with  II  and  X  increase  with  OJ.  The  cross 
products  then  automatically  have  the  same  signs. 

TOODY  has  no  built-in  provision  for  triangular  zones.  Tri¬ 
angular  zones  can  be  set  up  but  will  fail  the  time  step  calculations 
in  SETUP  1  and  MAIN  and  the  momentum  routine  in  MAIN.  If  triangular 
zones  are  to  be  used  then  these  coding  areas  must  be  changed.  The 
time  step  can  be  computed  using  SQRT(AL)  instead  of  distance.  The 
momentum  calculation  changes  required  depend  on  the  problem.  There 
is  no  general  change  which  will  work  for  all  possibilities. 

TOODY  cannot  be  successfully  run  (on  any  but  the  simplest  problems) 
by  a  casual  user  --  e.g.,  one  who  spends  a  day  or  two  reading  the 
documentation  and  then  inputs  his  problem.  The  code  has  not  been 
written  by  anyone  which  will  allow  this  type  of  black  box  treatment. 

To  consistently  run  this  code,  or  any  hydrocode,  successfully  requires 
patience  and  persistence  and  a  willingness  to  dig  into  the  detailed 
coding. 
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APPENDIX  A 


PLAIN  CONCRETE  CONSTITUTIVE  RELATIONS 
AT  HIGH  STRESS  LEVELS 

Calculations  of  projectile  impacts  or  explosive  detonations  in  concrete 
require  a  model  for  concrete  constitutive  relations  at  high  stress  levels. 

The  purpose  of  this  appendix  is  to  present  a  model  designed  to  simulate 
i  j.icrete  response  in  the  1Kb  to  300Kb  stress  region  and  be  reasonably  accu¬ 
rate  in  the  lower  stress  regions.  It  is  based  on  a  very  small  amount  of 
data.  Most  concrete  testing  has  been  conducted  at  stress  levels  below 
10,000  PSI  (less  than  1  Kilobar).  There  is  a  great  amount  of  research  still 
needed  to  adequately  answer  many  basic  questions  concerning  concrete  response 
at  higher  stress  levels. 

The  model  will  be  constructed  primarily  from  Hugoniot  data  on  one  specif¬ 
ic  concrete^  and  static  yield  strength  data  on  several  concretes^. 
Assumptions  will  be  made  for  a  Poisson's  ratio  value,  a  shear  modulus  formula, 
a  plastic  flow  rule  and  other  parameters. 

Gregson  generated  Hugoniot  dat~  to  stress  levels  over  500  Kilobars  for 
one  specific  concrete.  He  used  a  granite  aggregate  with  an  average  particle 
diameter  of  1/8  inch.  The  initial  density  varied  from  mix  to  mix,  but  aver¬ 
aged  2.185  gm/cc.  The  solid  density  had  an  average  value  of  2.673  gm/cc. 

The  concrete  consisted  of  18%  voids,  25%  feldspar  grains,  20  to  25%  quartz 


(1)  Gregson,  "A  Shock  Wave  Study  of  Fondu-Fyre  WA-1  and  a  Concrete  '.  i'N 
2797F,  Feb  72. 

(2)  Chinn  and  Zimmerman,  "Behavior  of  Plain  Concrete  Unde--  ,  ».  ■ 
Compression  Loading  Conditions",  AFWl  TR  64-163,  A ,r:  •• 
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grains  and  25  to  30%  cement  paste  by  volume.  The  small  aggregate  size  was 
dictated  by  his  shock  generating  method,  which  consisted  of  gas-  un  launched 
flyer  plates.  He  employed  several  measurement  techniques  to  establish 
stress  and  velocity.  His  final  data,  converted  to  stress  vs  excess  com¬ 
pression,  appears  in  Figures  A-l  and  A-2.  The  solid  line  is  simply  a  curve 

which  appears  to  fit  the  higher  stress  data  reasonably  well  and  also  takes 

*■ 

into  account  an  expected  quartz  phase  transition  at  150  Kb.  The  curve  for 
the  Hugoniot  is  drawn  through  the  upper  data  points  for  the  praqmatic  reason 
that  a  lower  curve  does  not  seem  to  adequately  predict  stress  levels  occur¬ 
ring  in  impact  tests. 

The  lowest  point  in  Gregson's  data  is  at  0.75  Kb.  This  is  a  fairly 
consistent  precursor  value  seen  in  his  stress  vs  time  data.  That  stress 
level  will  be  used  to  define  the  point  at  which  significant  cracking  and 
void  filling  occurs. 

The  unloading  curve  in  the  figures  is  simply  a  path  which  allows  the 
concrete  to  return  to  its  solid  density  when  loaded  past  the  point  where 
all  voids  are  filled.  It  is  a  straight  line  which  connects  with  the  loading 
curve  when  the  slopes  of  the  two  curves  are  equal.  There  is  no  data  for 
this  line  beyond  the  probable  solid  density  point  at  zero  stress. 

Measured  values  of  Gruneisen's  ratio  for  concrete  vary  but  are  typically 
around  0.1.  This  small  value  will  be  ignored  and  Gregson's  data  will  be  used 
to  predict  pressure  for  any  density  and  internal  energy  values. 

The  yield  strength  of  concrete  increases  as  confining  pressure  is  in¬ 
creased.  The  only  data  available  at  high  confining  pressures  is  static  data 
generated  by  Chinn.  Chinn's  results  are  shown  in  Figures  A-3  and  A-4  along 


FIGURE  A- 2  -  HUGONIOT  STRESS  DATA 


*ith  some  typical  low  stress  level  results^3).  In  the  figures  Y  is  defined 
as: 

Y  =  [((^  -  o2)2  +  (o2  -  a3)2  +  (a]  -  o3)2)/2]1/2 
Pressure  is  defined  as: 

P  =-1/3  (o"j  +  o2  +  03) 

whe>  the  0-j  are  principle  stresses  at  failure  and  are  measured  positive  in 
4  itsion.  This  value  of  pressure  differs  from  the  actual  confining  pressure 
generated  hydraulically  by  Chinn.  However,  it  is  consistent  with  the  way 
in  which  it  is  typically  defined  in  hydrocodes.  Values  of  Y  and  P  are 
normalized  in  the  figures  by  fc,  the  unconfined  compressive  strength  of  the 
concrete.  Chinn  used  mixes  with  fc  values  varying  from  3500  PSI  to 
10,000  PSI. 

Two  "theoretical"  points  are  shown  on  Figure  A-3.  One  describes 
failure  at  the  unconfined  compressive  strength  point— i.e.,  Y  *  fc  and 
P  =  1/3  f£.  The  other  is  included  because  a  minimum  pressure  criterion  will 
be  used  to  define  tensile  failure.  This  minimum  pressure  was  chosen  as  that 
occurring  in  an  unconfined  tensile  test  assuming  failure  when  0  =  0.1  fc  (*). 
It  is  numerically  pleasing  to  have  Y  =  0  at  the  point  where  P  =-0.1  fd/3. 

This  tensile  treatment  is  reasonable  for  a  model  not  expected  to  be  used 
extensively  at  low  stress  levels. 

The  yield  strength  curves  of  Figures  A-3  and  A-4  can  be  approximated  with 
three  straight  line  segments: 

(3)  McHenry  and  Kami,  "Strength  of  Concrete  Under  Combined  Tensile  and 
Compressive  Stress",  AC  I,  54,  10,  Apr  58. 

(4)  Crawford,  Higgins  and  Bultmann,  "The  Air  Force  Manual  for  Design  and 
Analysis  of  Hardened  Structures",  AFWL  TR  74-102,  Oct  74. 
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Y  =  3  (P  +  0.1  fc/3)/l .1  for  -0.1  fc/3  <  P  <  fc/3 

Y  *  P  +  2/3  fc  for  fc/3  <  P  <  30  fc 

Y  =  30.67  fc  for  P  >  30  fc 

The  limit  at  P  =  30  fc  is  not  indicated  in  Chinn's  data.  It  is  included  here 
because  it  is  felt  that  there  should  be  a  limiting  value  based  on  similar 
data  for  soils  and  rocks.  It  is  included  at  the  level  shown  simply  because 
it  is  known  to  not  exist  at  a  lower  level. 

It  should  be  emphasized  that  the  yield  strength  data  was  developed  in 
static  tests.  At  higher  strain  rates  the  yield  strength  may  increase  sub¬ 
stantially.  The  static  data  is  used  because  adequate  high  strain  rate 
data  is  not  available. 

Now  that  Y  has  been  determined  the  Hugoniot  stress  data  generated  by 
Gregson  can  be  reduced  to  pressure  data.  This  will  be  accomplished  by  in¬ 
voking  the  traditional  relationship 

a  =  P  +  2/3  Y 

which  hopefully  is  sufficiently  valid  for  uniaxial  tests  in  concrete. 

Letting  fc  =  5000  PSI  =  0.345  Kilobars,  we  can  use  the  yield  strength 
curve  to  define  Hugoiniot  pressure  points.  These  are  shown  for  loading  in 
Figure  A-5  by  the  dashed  line.  The  difference  between  the  stress  and  pres¬ 
sure  curves  builds  up  to  a  maximum  of  7.05  Kb  (2/3  x  30.67  x  fc)  at  a 
pressure  of  10.35  Kb.  The  unloading  curve  is  left  untouched- -remember  there 
is  no  data  but  the  solid  density  point  to  define  it  anyway. 

Figure  A-6  is  a  plot  of  the  resulting  P  -  y  curve  to  166  Kb.  It  should 
be  noted  that  there  is  really  no  data  between  63  and  166  kilobars.  The 
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FIGURE  A-5  -  HUGONIOT  STRESS  AND  PRESSURE 


y  =  p/po  -  1 


FIGURES  A- 6  -  HUGONIOT  PRESSURE 


phase  transition  is  consistent  with  data  at  and  above  the  166  kilobar 
level,  but  it  may  not  be  drawn  correctly  in  every  detail. 

Below  the  precursor  pressure  of  0.358  Kb  (stress  level  0.75  Kb), 
pressure  can  be  fitted  with  a  straight  line: 

P  =  KoP 

Kq  i  determined  from  the  Young's  Modulus  empirical  equation^). 

E  =  33  W  1 *5  (fc)®*^  PSI  (W  =  density  in  lb/ft^)  and  the  relationship 
K  =  E/3  (1  -  2  v) 

Poisson's  ratio,  v,  was  chosen  as  0.2. 

Designating  the  value  of  y  at  0.358  Kb  as  pE,  the  loading  pressure 
curve  can  be  fit  with  the  following  equations: 

P  (KB)  RANGE  OF  y 

KqP  •  -oo  to  PE 

0.358  +  78.62  (y-yE)  pE  to  0.1  • 

8  +  130  (p-0.1)  0.1  to  0.2 

21  +  420  (p-0.2)  0.2  to  0.3 

63  +  650  (p-0.3)  0.3  to  0.36607  (lockup  point) 

76  +  625.5  (p-0.32)  +  1720.19  (p-0.32)2  0.36607  to  0.414 

150  +  39.68  (p-0.414)  0.414  to  0.54  (transition  region) 

166  +  360.5  (y-0.572)  +  864  (y-0. 572)2  0.54  to  +  00 

The  maximum  density  unloading  curve  is  fit  by 
P  =  784  (p-0.223) 

For  pf/|  (the  largest  value  of  y  seen  by  the  concrete)  values  less  than  uE, 
loading,  unloading  and  reloading  take  place  along  the  same  path. 

For  pty  greater  than  the  value  at  lockup  (0.36607)  loading,  unloading  and 
reloading  take  place  along  the  same  path.  For  yu  values  in  between  these 


extremes  unloading  and  reloading  up  to  the  loading  curve  take  place  along 
Ky  lines  where  K  is  chosen  to  vary  linearly  between  Kq  and  784  Kilobars. 

These  loading  curves  and  some  unloading/reloading  paths  are  shown  in 
Figures  A-7  and  A-8. 

The  shear  modulus  6  is  calculated  assuming  that  the  Hooke's  law 
relationship  holds.  That  is, 

G  =  3M-2y)  _  ic 

2fr+v) 

where  Kk  is  the  bulk  modulus,  p  _  9.P  _. 

Deviators  which  cause  the  yield  surface  criterion  to  be  violated  are 
reset  in  the  same  manner  as  for  a  VonMises  surface.  It  can  be  shown  that 
the  resulting  plastic  flow  rule  is  inconsistent  with  a  Mohr-Coulomb  yield 
surface  such  as  proposed  for  concrete.  However,  the  alternative  is  a  very 
complex  computational  model.  The  extra  complexity  cannot  be  justified  con¬ 
sidering  the  source  of  the  yield  strength  model.  Strain  rate  effects  are 
probably  far  more  important  than  a  consistent  flow  rule. 

There  are  many  deficiencies  in  the  model.  However,  it  uses  the  avail¬ 
able  data  in  a  reasonable  manner.  Any  significant  shortcomings  in  the  model 
can  be  corrected  only  when  more  data  is  available  on  many  types  of  concretes 
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